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ABSTRACT 


This  research  project  was  undertaken  to  evaluate 
the  applicability  of  queuing  theory — in  particular  the 
priority  queuing  theory — to  the  evaluation  of  priority 
queuing  situations  in  military  information  systems. 

The  queuing  theory  has  been  applied  to  the  evalu¬ 
ation  of  the  queuing  situations  in  the  473-L  System. 
The  waiting-time  distribution  for  a  single-server, 
head - o f - th e -  1 i ne ,  priority  queuing  model  has  been 
evaluated.  The  application  of  “Variance  Reduction 
Method”  to  the  Simulation  of  Priority  Queuing  Systems 
has  been  investigated  Finally,  guides  and  procedures 
for  the  application  of  queuing  models  in  structuring 
information  systems  were  outlined. 
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I  INTRODUCTION 


A.  OBJECT  I VES 

The  objectives  of  the  research  effort  described  in  this  report  are: 

(1)  To  evaluate  the  applicability  of  queuing  theory  to 
evaluation  of  priority  queuing  situations  in  military 
information  systems. 

(2)  Develop  additions  and/or  modifications  to  the  ex¬ 
isting  priority  queuing  techniques  —  those  that  are 
deemed  necessary  from  the  application  studies — in 
sufficient  detail  for  use  as  a  system  design  aid. 

(3)  To  establish  guides,  procedures,  and  principles  for 
the  application  of  appropriate  queuing  models  in 
structuring  information  systems. 

B.  SCOPE 

The  473-L  Simulation  Model  has  been  chosen  as  the  sample  situation 
for  use  in  the  evaluation  of  the  applicability  of  queuing  theory  to  the 
evaluation  of  priority  queuing  situations  in  military  information  systems. 
The  analytical  analysis  of  the  473-L  Model  is  presented  in  Appendix  A. 
Results  of  the  analytical  investigation  have  been  compared  with  that  ob¬ 
tained  from  simulation  and  are  presented  in  Sec.  II. 

Usually  the  analysis  of  priority  queuing  is  limited  to  the  mean  of 
the  operational  measures.  An  analytical  expression  for  the  priority 
queuing  waiting-time  distribution  has  been  developed  during  the  previous 
research  effort.  It  is  shown  in  Sec.  Ill  that  this  waiting-time  distri¬ 
bution  can  be  approximated  by  a  simpler  expression  when  the  system  loading 
and  the  relative  loading  of  the  higher  priority  to  the  lower  priority 
customers  meet  certain  requirements. 

Although  it  is  possible  to  analytically  investigate  the  queuing 
situations  in  information  systems,  this  usually  involves  extensive  ideali¬ 
zations  of  the  system  parameters.  The  results  obtained  are  usually 
limited  to  the  mean  of  the  measures.  This  has  been  shown  in  Sec.  II.  For 
more  complex  information  systems,  where 
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(1)  the  rnter -arrival  interval  distributions  or  the  service 
time  distributions  are  not  exponential,  and 

(2)  where  queuing  discipline  is  of  the  priority  type  and  the 
service  center  contains  multiple  parallel  servers,  the 
only  practical  technique  for  system  evaluation  is  the 
simulation  modeling  technique. 

It  is  therefore  highly  desirable  to  develop  simulation  strategies  that 
can  be  used  to  increase  the  effectiveness  of  this  technique.  Three  as¬ 
pects  of  the  simulation  strategies  have  been  investigated  in  Sec.  IV. 
These  are  the  stopping  rules  and  the  chopping  rules  of  the  simulation 
experiment,  and  the  application  of  “Variance  Reduction  Methods9’  to  the 
Simulation  of  Priority  Queuing  Systems. 

finally,  the  guides  and  procedures  for  the  application  of  the 
Queuing  Models  are  outlined  in  Sec.  V. 
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II  RESULTS  01  ANALYTICAL  ANALYSIS  OF  47 3 -L  SYSTEM 


A.  INTRODUCTION 

The  principal  purpose  of  the  queuing  analysis  of  the  473-L  system 
is  to  investigate  the  possibility  of  applying  the  queuing  theory  to 
obtain  some  of  the  queuing  measures  of  a  real  information  system.  The 
473-L  system  model  has  been  chosen  for  this  purpose  because  a  simulation 
model  of  this  system  has  been  constructed  and  simulated  results  of  this 
system  exist.1  Thus,  it  is  possible  to  check  the  analytical  results  with 
the  simulated  results. 

B.  DESCRIPTION  OF  473-L  SYSTEM 

In  this  section  a  brief  description  of  the  473-L  Queuing  model  is 
given.  For  a  more  detailed  description  see  Appendix  A  or  Reference  1. 
Figure  III  shows  the  473-L  Queuing  Simulation  Model.  Requests  for 
service  arrive  at  Buffer  1,  Q- 1 ,  from  the  source  of  data.  The  requests 
are  categorized  into  three  priority  classes.  The  queuing  discipline  at 
Queue  1  is  that  of  interrupt  priority  type.  All  requests  require  an 
initial  service  by  the  computer.  Some  of  the  requests  require  further 
service(s)  from  the  computer.  Each  further  service  by  the  computer  is 
preceeded  by  a  disc  service  which  is  to  retrieve  the  desired  information. 
Those  requests  that  require  further  service(s)  from  the  computer  would 
have  to  enter  the  Buffer  2,  Q- 2,  to  wait  for  the  retrieval  operation. 

The  queuing  discipline  at  Q-2  is  of  the  ordered  queue  type.  When  the 
required  number  of  computer  services  have  been  completed,  an  output 
message  is  sent  to  the  display. 


Reference*  *re  lifted  *t  the  end  of  the  report. 
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FIG.  11-1  473-L  SIMULATION  MODEL 
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C.  RESULTS  OF  THE  ANALYSIS 


The  analytical  expressions  for  some  of  the  queuing  measures  are 


(1) 

the 

utility  of  the  computer  pc 

(2) 

the 

utility  of  the  disc,  pd 

(3) 

the 

mean  response  time  of  all 

requests,  Tr 

(4) 

the  mean  response  time  of  each 
requests ,  T  ,  and 

priority  class  of 

(5) 

the 

mean  storage  size  for  Q- 1 

and  Q- 2 ,  Ly 1  and  L 

The  utility  of  the  computer  is  defined  as  the  percentage  of  time 
the  computer  spent  in  servicing  requests.  Similarly,  the  utility  of  the 
disc  is  the  percentage  of  time  that  the  disc  spent  in  servicing  requests. 
The  mean  response  time  is  defined  as  the  mean  elapsed  time  between  arrival 
of  the  new  request  at  Q-  1  to  the  time  when  it  leaves  the  computer  after 
final  service  by  the  computer.  The  mean  storage  size  is  defined  as  the 
mean  number  of  words  in  the  buffers. 

Detailed  derivations  of  these  expressions  are  contained  in 
Appendix  A.  The  expression  for  each  of  the  above  queuing  measures  are 
listed  below; 

(1)  The  utility  of  the  computer 
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(2)  The  utility  of  the  disc 
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(3)  The  mean  response  time  (sec)  of  all  requests 
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(4)  The  mean  response  lime  (sec)  of  l  lie  itli  priority  request 
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(5)  The  mean  storage  size  (words)  has  been  developed  under 
two  sets  of  storage  rule,  Buie  1  and  Buie  2.  See 
Appendix  A  for  more  detailed  explanation  of  these 
rules. 

The  mean  storage  size  for  Q- 1  and  Q- 2,  aeeording  to  Buie  1,  are 
given  as 
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The  mean  storage  size  for  Q- 1  and  Q- 2,  aeeording  to  Buie  2,  are 
given  as 
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a  nd 
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The  mean  arrival  rate  of  requests  from  the  source 
of  data. 

The  mean  service  time  of  the  computer. 

The  percentage  of  new  requests  that  require  no 
disc  search. 

The  housekeeping  time  due  to  an  interrupt. 

The  probability  that  a  request,  upon  its  arrival 
at  Q-  1 ,  finds  that  the  request  that  is  being 
serviced  by  the  computer  has  a  lower  or  equal 
priority. 

The  mean  service  time  for  the  disc  search 
request. 

The  mean  service  time  for  the  data  disc  search 
request. 

The  mean  service  time  for  the  program  disc 
search  request. 

The  mean  number  of  disc  searches  required  by  a 
new  request. 

The  percentage  of  new  requests  that  are  ith 
priority  requests. 

The  mean  length  of  a  new  request. 

The  mean  length  of  a  disc  search  address. 

The  mean  length  of  the  message  retrieved  by  a 
data  disc  search  request. 

The  mean  length  of  the  message  retrieved  by  a 
program  disc  search  request. 

The  proportion  of  n  disc  searches  that  are  the 
data  disc  searches. 

The  mean  length  of  a  display  message. 


The  results  of  the  simulation  for  twenty-four  cases  are  listed  in 
Table  I,  as  extracted  from  Reference  1. 
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The  queuing  measures,  computer  utility  Pc,  the  disc  utility  pd,  and 
the  response  time  T  ,  are  computed,  using  Eqs.  (1),  (2),  and  (3),  respec¬ 
tively,  for  the  twenty-four  cases.  The  computed  values  are  listed  in 
Tabic  11  together  with  the  simulated  values.  The  p  is  computed  on  the 
basis  of  simulated  mean  number  of  disc  search.  The  p ^  is  computed  on  the 
basis  of  simulated  computer  utilization  and  simulated  mean  number  of  disc 
search.  Finally,  the  simulated  n,  p c ,  and  pd  were  used  to  compute  the  T 
The  reason  for  using  the  simulated  n,  Pc ,  and  pd  in  computing  each  of  the 
queuing  measures  is  to  make  the  system  environments  of  the  analytical  model 
match  as  closely  as  possible  that  of  the  simulation  model. 

Table  III  is  a  listing  of  the  simulated  and  computed  mean  response 
time  for  the  ith  priority  request,  T . ,  for  the  above  cited  twenty-four 
cases.  Equation  (4)  was  used  to  compute  the  7V.  Once  again,  the  simu¬ 
lated  n,  pc  and  pd  were  used  in  the  computation. 

Equations  (5)  and  (6)  were  used  to  compute  the  mean  storage  size  in 
Q- 1  and  Q- 2  for  the  twenty-four  cases.  These  are  listed  in  the  second  and 
third  column  of  Table  IV.  The  fourth  column  is  the  sum  total  of  the  first 
two  columns.  The  Lol  and  Lq2  listed  in  the  fifth  and  sixth  columns  were 
computed  according  to  Eqs.  (7)  and  (8).  Their  sums  are  listed  in  the 
seventh  column.  The  simulated  results  are  listed  in  the  eighth  column. 

The  percent  deviation  of  the  simulated  results  from  the  computed  results, 
according  to  Rule  1  and  Rule  2,  are  tabulated  in  the  ninth  and  tenth 
col umns  . 

D.  DISCUSSION 

It  is  noted  in  Table  II  that  the  analytical  results  are  generally  in 
agreement  with  the  simulated  results,  except  for  a  few  specific  cases. 

The  most  noticeable  case  is  Case  18,  the  only  one  where  the  computed  pc  is 
higher  than  the  simulated  pc .  Upon  further  examination,  it  is  seen  Case  18 
is  the  only  one  where  the  disc  is  saturated.  The  computed  disc  utility  is 
142.36%,  which  should  be  interpreted  as  the  offered  load  at  Q- 2,  but  the 
actual  utility  of  the  disc  is  100%.  Since  the  disc  is  saturated,  the  mean 
rate  at  which  the  request  returns  to  the  Q- 1  after  the  disc  search  becomes 
independent  of  the  arrival  rate  of  new  requests.  This  rate  is  the  mean 
disc  search  rate,  which  is  approximately  1 / h d .  Thus,  the  mean  arrival  rate 
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CASE 

NO. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 


Table  III 


COMPUTED  AND  SIMULATED  MEAN  RESPONSE  TIME 


MEAN  RESPONSE  TIME  («ec) 

PRIORITY 

1st 

2nd 

lid 

Computed 

Simulated 

Coapu  ted 

Simulated 

Co  P  P  u  1 

Led  Simulated 

2.56 

2. 

25 

2.58 

2.39 

2 

.64 

2.54 

2.56 

2. 

25 

2.58 

2.39 

2 

.64 

2.54 

10.11 

11. 

29 

10.21 

10.33 

10 

.43 

10.38 

10.11 

11. 

29 

10.21 

10.33 

10, 

.43 

10.38 

2.34 

2. 

16 

2.40 

2.25 

2 

.56 

2.68 

15.  15 

12. 

26 

15.18 

12.28 

15, 

.29 

13.82 

11.  10 

12. 

51 

11.13 

11.59 

11 

.34 

11.82 

7.06 

6. 

10 

7.03 

6.93 

7, 

.11 

7.94 

22.96 

20. 

79 

23.03 

21.53 

23 

.13 

22.36 

16.  10 

15. 

59 

16.11 

15.78 

16, 

.  16 

17.81 

7.27 

5. 

86 

7.37 

6.49 

7, 

.69 

7.06 

3.01 

2. 

88 

3.08 

3.51 

3, 

.19 

3.33 

9.20 

7. 

01 

9.30 

9.05 

9, 

.89 

10.38 

4.39 

3. 

81 

4.84 

4.08 

7. 

,06 

7.76 

9.  11 

8. 

32 

9.33 

8.23 

10. 

.09 

9.48 

3.88 

3. 

55 

4.24 

3.87 

5. 

,75 

5.54 

4.47 

3. 

41 

5.00 

4.48 

8. 

,23 

8.00 

6066.00 

687. 

30 

6082.00 

692.50 

6119. 

.00 

646.70 

57.30 

42. 

60 

57.48 

51.80 

57. 

,83 

42.60 

8.56 

8. 

85 

8.58 

7.99 

8. 

63 

8.94 

13.64 

15. 

97 

13.66 

14.40 

13. 

,72 

14.56 

2.86 

2. 

88 

2.95 

3.04 

3. 

,23 

3.33 

7.86 

8. 

41 

8.02 

6.94 

8. 

56 

8.37 

7,04 

5. 

54 

7.08 

7.28 

7. 

19 

8.00 

11 


Table  IV 


(DM  1*1  TED  AM)  SIMULATED  MEAN  (X)RE  LENGTH 


CASE 

NO. 

MEAN  CORE  LENGTH  (WORDS) 

PERCENT 

DEVIATION 

100  U  -  L  )/L 

Compu  te d 

Sim  u 1  a  t  ed 

L 

s 

Rule  1 

Rule  2 

According  to 

Rule  1 

According  to  Rule  2 

lq\ 

LQ2 

Total  — L 

c 

lqi 

LQ2 

Total  — L 

1 

460 

680 

1,  140 

465 

695 

1,160 

651 

42.9 

43.9 

o 

1030 

1.610 

2,640 

1035 

1  625 

2,660 

1,700 

35.6 

36.1 

3 

520 

1,790 

2,310 

535 

1,845 

2,  380 

3,400 

-47.2 

-42.7 

4 

1170 

4,110 

5,280 

1180 

4,170 

5,350 

7,500 

-42.0 

-40.2 

5 

1160 

1,470 

2,630 

1170 

1,475 

3,645 

1,900 

27.8 

28.2 

6 

1010 

3,370 

4,  380 

1020 

3,370 

4,390 

5, 100 

-16.4 

-16.3 

7 

990 

2,360 

3,350 

995 

2,385 

3,380 

3,700 

-10.4 

-  9.4 

8 

950 

1  460 

2,410 

960 

1,465 

2,425 

1,400 

42.0 

42.2 

9 

970 

3,370 

4,  340 

980 

3,380 

4,360 

5,800 

-33.7 

-33.0 

10 

96  0 

2,330 

3,  290 

965 

2,335 

3,  300 

3,500 

-  6.4 

-  6.  1 

11 

545 

1,250 

1,795 

565 

1,300 

1,865 

2,000 

-11.4 

-  7.3 

12 

505 

715 

1,220 

520 

7  40 

1,260 

1,000 

18.0 

20.5 

13 

595 

1,  530 

2,  125 

620 

1,605 

2,225 

3,000 

-39.3 

-34.9 

14 

1010 

745 

1,755 

1060 

810 

1,870 

1,700 

3.  1 

9.0 

15 

670 

1,  520 

2, 190 

700 

1,610 

2,310 

2,900 

-32.4 

-25.4 

16 

735 

665 

1, 400 

780 

710 

1,490 

1,200 

14.3 

19.5 

17 

965 

690 

1,655 

1040 

755 

1,795 

1,700 

-  2.7 

5.3 

18 

1060 

1 , 342,000 

1,343,000 

1075 

1,356,000 

1,357, 100 

613,000 

52.9 

54.8 

19 

1050 

12,880 

13,930 

1060 

13,010 

14,070 

25,000 

-79.5 

-77.7 

20 

970 

1,805 

2,775 

975 

1,815 

2,790 

2,200 

20.7 

21.2 

21 

950 

1,960 

2,910 

955 

1,965 

2,  920 

2,600 

10.7 

11.1 

22 

1320 

1,930 

3,250 

1335 

1,970 

3,305 

3,500 

-  7.7 

-  5.9 

23 

605 

1,315 

1,920 

635 

1,  380 

2,015 

2,400 

-77.1 

-68.8 

24 

480 

1,295 

1,775 

495 

1,330 

1,825 

2,300 

-29.6 

-26.3 

of  new  requests  and  return  requests  is  \  +  l//i^.  It  follows  that  the  com* 
puter  utility  due  to  the  analysis  of  requests  is: 


P\  '  <*■  +  1 /hd)h( 


(9) 


As  may  be  seen  from  Eq 
Pc  ~  1,  which  is  the  worst 
to  interruption  p**  and  the 
requests  p '  is: 


s.  (A- 3)  and  (A -4),  under  the  assumption  of 
case,  the  ratio  of  the  computer  utility  due 
computer  utilization  due  to  the  analysis  of 


p"c  hk[  1  ♦  ( N  -  1  )Z] 

K  ’  M 


(10) 
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For  tire  system  under  study,  X  =  1/4,  Z  =  0.76,  and  hh  *  1  msec;  for 
Case  18  where  hc  *  50  msec,  the  ratio  p"c/ p'c  =  0.015,  which  means  Pc  is 
small  in  comparison  to  p\  Thus,  Pc  may  be  approximated  by  Pc.  More 
specifically,  under  disc  saturation  conditions,  pc  may  be  approximated  by 
Eq.  (9).  For  Case  18,  where  X  -  1/3  requests/sec  and  hd  *  0.36  sec/request, 
p'c  =  pc  =  15.55%,  which  compares  favorably  with  the  simulated  result  of 
15.31%.  Thus,  under  disc  saturation  conditions  Eq .  (9)  instead  of  Eq .  (1) 

should  be  used  to  compute  the  utility  of  the  computer. 

The  computed  Tr  for  Case  18  is  about  nine  times  the  simulated  Tr ;  for 
Case  19,  the  computed  Tr  is  about  16%  higher  than  the  simulated  Tf.  In 
both  of  these  cases,  the  disc  utility  is  near  or  at  100%.  It  is  assumed 
the  simulation  started  from  the  empty  state,  with  no  request  in  the  system, 
and  then  proceeded  to  build  up  to  its  steady  state  level.  In  the  case  of 
100%  utility  there  is  no  steady  state  level  since  the  queue  at  the  disc 
will  continue  to  build  up  without  bound.  The  time  of  buildup  could  be 
rather  long  when  the  utility  of  the  server  is  near  100%.  The  simulated 
results  are  based  on  1000  reply  messages  being  generated  within  the  system 
and  sent  out;  therefore,  it  is  conjectured  that  the  simulated  Tr  for 
Cases  18  and  19  does  not  represent  the  steady  state  T  .  They  are  heavily 
influenced,  especially  the  Tr  of  Case  18,  by  the  initial  condition  of  the 
system.  It  is  concluded  that  either  a  longer  simulation  run  is  needed 
for  situations  like  Cases  18  and  19,  or  that  means  should  be  devised  to 
eliminate  the  bias  caused  by  the  transient  condition. 

The  computed  Pc  are  less  than  the  simulated  Pc  in  all  cases,  with  the 
single  exception  of  Case  18,  as  was  noted  above.  One  possible  explanation 
is  the  bias  introduced  in  the  simulated  mean  arrival  rate  and  service 
time,  as  evident  in  the  bias  introduced  in  the  mean  disc  search.  In  all 
cases,  the  simulated  mean  disc  search  is  less  than  the  specified  mean  disc 
search. 

The  computed  mean  response  times  for  each  of  the  three  priority 
classes  seems  to  be  in  general  agreement  with  the  simulated  results  (see 
Table  III),  except  in  the  Cases  18  and  19.  The  above  remarks  on  these 
two  cases  apply  here  too. 

It  is  to  be  noted  in  Cases  3,  4,  7,  19,  20,  21,  and  23  that  the 
simulated  mean  response  time  for  the  higher  priority  messages  is  greater 
than  that  of  the  lower  priority  messages.  This  inconsistency  is  probably 
due  to  sampling  errors  of  the  simulation  process. 


13 


Tin*  simulated  mean  tort*  length  I'm  t*«it  h  case,  ns  shown  in  Table  IV, 
it*  p  r  (‘st*  n  t  s  a  single  sample  observation.  1‘lie  re  Tore,  it  is  rather  dill  mil  t 
to  indicate  the  degree  of  agreement  between  the  analytical  models  and  the 
simulated  model.  However,  an  observation  may  be  made  on  the  percentage 
deviation  of  the  simulated  mean  from  the  analytical  mean.  Table  V  shows 
the  proportion  of  cases  that  deviate  less  than  x  percent  for  Storage 
Hale  1  and  2.  Also  shown  is  the  normal  distribution  with  cr  =  1/3  mean. 


Table  V 

CUMULATIVE  DISTRIBUTION  OF  PERCENT  DEVIATION 
OF  SIMULATED  MEAN  CORE  LENGTH 
FROM  THE  ANALYTICAL  MEAN 


PERCENT 

PROPORTION  OF  CASES 

NORMAL 

DEVIATION 

Rule  1 

Rule  2 

DISTRIBUTION 

1  10 

16.7 

25.0 

23.6 

1  20 

41.7 

41.7 

45.  1 

^  30 

54.  1 

62.5 

63.2 

1  40 

70.9 

83.4 

77.0 

1  50 

87.5 

91.7 

8b .  6 

1  oo 

91.7 

96.0 

92.8 

i  70 

91.7 

9b. 0 

96.4 

1  80 

100.0 

100.0 

98.4 

A  comparison  of  the  observed  percent  deviation  distribution  with  the 
normal  distribution  seems  to  indicate  the  simulated  means  have  a  normal 
distribution  in  relation  to  the  analytical  mean.  This  implies  a  reason¬ 
able  agreement  exists  between  the  analytical  means  and  the  simulated  means. 

E .  CONCLUSIONS 

The  analytical  analysis  of  the  473-L  System  indicates  it  is  possible 
to  apply  elementary  queuing  theory  to  analytically  model  a  relatively 
complex  information  system,  and  to  obtain  some  of  the  queuing  measures  of 
the  svstem.  The  validity  of  the  analytical  model  may  be  verified  by  sim¬ 
ulation.  Once  verified,  the  analytical  model  may  be  used  by  the  system 
designer  in  considering  the  trade-off  of  system  elements. 
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Ill  NUMERICAL  EVALUATION  OF  WAITING-TIME  DISTRIBUTION 


A.  INTRODUCTION 

As  seen  in  Sec.  II,  by  suitable  idealization  of  the  information  sys¬ 
tem  we  can  apply  the  queuing  theory  to  obtain  the  means  of  the  queuing 
measures  of  the  system.  The  information  system  designer  is  often  inter¬ 
ested  in  the  distributions  of  the  queuing  measures  so  that  end  point  con¬ 
ditions  of  the  system  can  be  analyzed.  However,  analytical  expressions 
for  the  distributions  of  priority  queuing  measures  usually  are  not  avail¬ 
able,  or  if  available  they  seldom  are  in  readily  usable  form.  Thus,  an 
approximate  analytical  expression  that  is  readily  computable  is  a  welcome 
addition  to  the  information  system  design. 

During  the  previous  research  effort,  Eq.  (8)  of  Reference  2,  the 
waiting-time  distribution  of  the  second - pri ori ty  customer,  W 2(T) ,  in  an 
exponential  head-o f - the  - 1 i ne  priority  queuing  system  has  been  shown  to  be 


W2{>  T) 


1  -  P 

2v 


/( r )dr , 


for  p2  >  px 

(11a) 


1  -  P 

2v 


r  2 

/  / ( r  )dr , 


for  p2  <  /)j 


(lib) 


where 

p  =  the  loading  factor  of  the  system 

p  j  =  the  loading  factor  of  the  system  due  to  the  first 
priority  customers 

T  *  unit  of  time  in  units  of  mean  service  time,  p. 


a 


l 


a 
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a 


(p  -/>,)(  1  -  p) 


r 


fir) 


-  rT 


/(a2  -  r)(r  -  a,) 


[r(r  -  a)] 


Cox'  has  derived  a  simple  expression  for  the  distribution,  namely 


»2  ( >  T) 


-  (  \-P  .  )  (  \-P  )  T 

=  pe 


(12) 


But  Cox’s  expression  is  not  valid  for  all  ranges  of  p  j  value.  Since 

Eq .  (12)  is  a  relatively  simple  equation  it  is  of  interest  to  investigate 

II 

the  range  of  p^  value  where  W ^(>  T)  is  valid. 

It  is  seen  that  when  p j  <<  p,  then  ( p 2  -  p ^) / p  -  pj)  p,<X  —  1  -  p 
and  f  2  f(r)dr  —  0.  Therefore,  Eq.  (11)  becomes  that  of  Eq.  (12).  It  is 

also  to  be  noted  that  when  p  <<  1  implies  p.  <<  1,  which  means 

f(r)dr  —  0.  Thus,  Eq .  (11)  again  tends  toward  Eq.  (12). 

B.  RESULTS 

A  fttiffterical  evaluation  of  Eq.  (11)  was  carried  out  for  p  *  0.1,  0.2, 

OJ,  0,5,  0,7,  0,05,  a  fid  0,95,  and  for  pjp  *  0.1,  0.2,  0.5,  and  0.9. 

Ill*4  result**  Of  thfe  evaluation  art*  plotted  in  Fig.  I1I-1  through  6  for 

each  value  of  cxre.pt  -  U.J.  When  /•  r  0.1  Eq .  (11)  a. id  Eq .  (12)  are 
almost  identical.  The  abM  »  s.t  of  t  lie  figures  is  in  units  of 
V  ^20(1  -  p ) ( 1  -  p  j ) T l  ’ ( / n  f  -  In  L  0  ~  M  .  Th is  unit  of  A  i s  r hose  n  so 

that  Lq  .  (12)  appears  js  ?  ingle  n  r  v  c  for  a  specific:  value  of  /  ■  ^  .  In  the 

figure^  t  h»-  unmarked  straigh*  line  is  t  lie  wu  i  \  i  ug  -  l  i  me  distribution  due 
t-.  Cos.  If  ,  T*  . 

We  ♦  i'  ft  ni(i  the  fig  u  I  e  s  t  )i  a  t  t  h  <■  vh  l  t.  i  ng-  .  l  me  i  i  :  t  i  bu  t  i  on  due  t  O 

( .<■  \  approximates  that  of  Eq  .  \ll)  quite  well  ior  ,  }  ^  0 . 2p  and  for  a  1  J 

'  •  '  1  1  t  -  1  1  •  e\c  ept  i  on  when  p  i  s  nea  r*  0 . 5  .  When  f  is 

n  j 

near  0.5,  the  percentage  deviation  ol  ffgl^Dlrom  ft  2  (>  T)  becomes  rela¬ 
tively  large  at  the  extreme  value  of  Xt  as  much  as  45%  at  X  «  20.  As  the 
ratio  of  Pj/p  Increases,  the  percentage  deviation  becomes  larger,  which 
verified  the  analytical  observation  made  above. 
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Krom  the  figures  we  note  that  tq .  (12)  yields  a 

estimate  of  waiting -time  at  the  lower  value  of  X  and 
estimate  at  higher  value  of  X.  The  cross-over  point, 

i  n 

W 0 ( >  T )  =  H  0 ( >  T ) ,  varies  with  the  value  of  p . 


more  optimistic 
a  more  pessimistic 
where 
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W(X)  =  PROBABILITY  OF  EXCEEDING  "x"  DELAY  UNITS  BY  THE  SECOND  PRIORITY  CUSTOMERS 


X  --  [20  (l-p)(l-p,)M*]/('n/»-ln  I0's) 

T  •  -  8  I  •  7  -  9 


FIG.  Ill-l  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY  —  p  »  0.2 
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W(X)  =  PROBABILITY  OF  EXCEEDING  "x"  DELAY  UNITS  BY  THE  SECOND  PRIORITY  CUSTOMERS 


0  2  4  6  8  10  12  14  16  16  20 


X  *  [20  (l-p)('-p1)M,]/('nP  "In  'O'3) 

T« -Sli 7 -4 

FIG.  111-2  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY— p  -  0.3 
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W(X)  S  PROBABILITY  OF  EXCEEDING  “x“  DELAY  UNITS  BY  THE  SECOND  PRIORITY  CUSTOMERS 


X  =[20  (1  -  p)  (t -pj  M»]/(ln/>- In  10  3) 

T9-  91  •  7  ~  6 


FIG.  111-3  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY  —  p  =  0.5 
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W(X)  =  PROBABILITY  OF  EXCEEDING  "X  "  DELAY  UNITS  BY  THE  SECOND  PRIORITY  CUSTOMERS 


X  1  [20  (l ~p)  (l-p,)M»]/(liV>-ln  IO'J) 

Ti-Jiir*  • 


FIG.  II 1-4  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY  —  p  =  0.7 
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W(X)  r  PROBABILITY  OF  EXCEEDING  "x"  DELAY  UNITS  BY  THE  SECOND  PRIORITY  CUSTOMERS 


TI-JIIT-7 


FIG.  MI-5  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY  —  p  =  0.85 


oo 


W[X)  s  probability  of  exceeding  "x"  oelay  units  by  the  second  priority  customers 


FIG.  111-6  WAITING-TIME  DISTRIBUTION  OF  THE  SECOND  PRIORITY  —  p  -  0.95 
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IV  SI MULAT ION  STRATEGY 


A.  INTRODUCTION 

In  many  information  systems  the  customers  (messages)  are  categorized 
into  priority  classes  and  their  order  of  service  is  in  accordance  with  a 
given  priority  queuing  discipline.  An  example  of  this  type  of  information 
system  is  the  473-L  system,  as  described  in  Sec.  II.  It  is  sometimes  pos¬ 
sible  to  rdealize  the  parameters  of  the  system  and  analyze  its  queuing 
measures  analytically,  as  was  done  in  Sec.  II.  Because  there  is  a  rela¬ 
tive  scarcity  of  analytical  results  on  the  priority  queue,2  however,  it 
is  not  always  possible  to  analyze  the  queuing  measures  without  extreme 
idealizations.  Such  extreme  idealizations  may  make  the  results  of  the 
analysis  unrealistic.  Under  such  a  circumstance  the  method  of  simulation 
frequently  is  used  to  analyze  the  queuing  measures  or  to  check  the  valid¬ 
ity  of  idealizations.  It  would  seem  advisable,  therefore,  to  develop  a 
guide  for  efficient  design  of  the  simulation  experiment. 

The  present  study  investigated  the  stopping  rules  and  the  chopping 
rules  of  the  simulation  experiment,  and  the  application  of  “Variance  Re¬ 
duction  Methods”  to  the  Simulation  of  Priority  Queuing  Systems. 

B.  STOPPING  RULES 

In  simulation  the  question  often  asked  is  how-  long  should  be  the 
simulation  experiment?  Alternately  the  question  is  when  to  stop  the  sim¬ 
ulation  experiment?  One  possible  answer  is  to  stop  the  simulation  exper¬ 
iment  when  the  observed  variation  in  the  estimates  reaches  certain  specified 
confidence  interval  sizes.  Based  upon  a  concept  of  Conway,4  the  confi¬ 
dence  interval  about  X  has  been  shown  to  be  (Appendix  B)  : 

X  ±  aS(X  ) 

n 

where  a  =  the  number  of  standard  deviations  in  an  (1  -  CL )%  confidence 
interval . 
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s-(X) 


S3iX)  1  ♦  P. 


n  1  -  p. 


sHx) 


rht?'  ^-"(7.)2] 


1  z  * 

n  i=l  ' 


1  £ 

X ,  =  -  Z  y 

*  L  j=  i  " 


«> 


the  j t h  observation  of  the 
queuing  measure  in  the  tth 
block  of  the  simulation 
experiment . 


(13) 


n  +  2  1  +  p 

(n-l)2l~  P 


p  =  C/[(n  -  1)S2(*)] 


"i1  (xtxl+1)  -  (n  -  nanl)un2) 

r  *  1 


n  1 


-  I  xt 

n  -  1  t-  l 


1  n-l 

x  0  *  -  s 

n2  n  -  1  t-1  f 


(14) 


Given  L  and  a,  it  is  possible  to  determine  the  confidence  interval 
about  X  n  at  any  value  of  n ,  where  n  is  the  length  of  simulation.  Since 
the  confidence  interval  can  be  computed  as  the  simulation  progresses,  it 
is  possible  to  stop  the  simulation  whenever  the  confidence  interval  is 
less  than  certain  specified  size. 
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Aliy  stopping  rnli1  which  depends  on  the  obs  e  r  v  n  i  i  oiis  t  Items  elves  can 
produce  bins.  For  example,  suppose  that  the  sample  mean  X ,  and  sample 
variance  5“(Ar)  are  i  ndepe  nden  t  1  y  distributed.  Tins  is  true  i  f  X  (  is 
normally  distributed.  Suppose  the  stopping  rule  is; 

Stop  when  *S(  An)  £  KX n 

Then  when,  fora  given  vS(A£),  Xn  is  low  the  simulation  will  be  stopped. 

If  X  is  high,  further  samples  will  be  taken  which  will  tend  to  reduce 
X  .  Thus  this  rule  favors  lower  values  of  X  and  is  biased. 

n  n 

For  X  t  normally  distributed,  the  rule: 

Stop  when  5(A£)  £  A 

is  unbiased  because  now  the  stopping  point  is  independent  of  X n . 

Another  example  is  X  t  with  an  exponential  distribution.  Here  S(X n) 
tends  to  be  proportional  to  X n .  »S  ( A'  n  )  =  (XAf  n  .  Th en  the  rule: 

Stop  when  *S  ( X  n  )  £  K 

is  biased  because  £  (A  )  <  K  implies  X  <  A  Ct  and  again  low  values  of  X 

n  —  *  n  - —  0  n 

will  satisfy  the  rule  and  stop  simulation.  In  this  case  the  rule: 

Stop  when  £(A'  )  <  KX 

1  n  —  n 


will  be  approximately  unbiased. 

The  X  t  can  usually  be  assumed  to  be  normally  distributed,  therefore 
the  unbiased  stopping  rule  should  be: 


Stop  when  S(Xn)  =  A 


C.  CHOPPING  RULES 

systems  is  the  selection  of 
selected  at  random  from  the 
is  in  steady-state  the  choice 
average  statistics.  A 


A  problem  in  the  simulation  of  queuing 
initial  state.  If  the  initial  state  is  not 
true  distribution  of  states  when  the  system 
will  introduce  some  bias  into  the  resultant 
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possible  course  of  action  to  avoid  bias  is  to  start  the  simulation  from 
empty  state  and  chop-off  the  first  part  of  the  run  from  the  average 
statist i c  s . 

A  chopping  rule  was  developed  and  is  presented  in  Appendix  C.  This 
chopping  rule  is  based  upon  the  assumption  that  the  simulation  experiment 
would  start  from  the  empty  state  and  progress  to  the  s teady - s ta te .  Hence 
the  parameter,  which  in  this  case  would  be  the  mean  waiting-time  as  esti¬ 
mated  from  the  initial  part  of  the  simulation,  would  be  lower  than  that 
obtained  from  the  later  part  of  the  simulation.  A  statistical  test  pro¬ 
cedure  that  can  test  the  hypothesis  that  the  parameter  as  estimated  from 
the  initial  part  of  the  simulation  is  less  than  that  from  the  later  part 
of  the  simulation  can  be  used  to  determine  the  point  of  chopping.  The 
testing  procedure  is  as  follows  (using  the  notations  of  the  stopping  rule) 

Let  the  sequence  of  the  estimates  of  the  parameter  from  the  simula¬ 
tion  be 


-  1  " 

.V  =  -  2  X,,  m  •  1,  2,  3 . n 

*  m  i-4  ' 


Starting  from  m  9  1  compute  the  statistics 


t 

m 


X 


n 


-  m 


-  X 


m 


l)S*Wa)  +  (n  -  m 
n  -  2 


1)S2(X 


-1  m 


(15) 


where 


X 


nX  -  mX 

n  m 


n  -  n 


S2(X 


n  -  m  -  1 


t(n  -  ns*(jr. )  -  (m  -  1  )S2(X  )  -  (n  -  m)X2  +  nX2  -  mX2] 


Find  m'  -  max  {m:tm  2)and  chop  at  . 

Based  upon  experience  one  may  ignore  the  chopping  rule  unless  the 
system  loading  factor  is  greater  than  0.9. 
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I).  Tilt:  APPLICATION  OF  "VARIANCE  REDUCTION  METHODS" 
TO  THE  SIMULATION  01  PRIORITY  QUEUING  SYSTEMS 


l.  General 

a .  I NTRODUCT I  ON 

The  history  of  variance  reduction  is  very  closely  associated 
with  that  of  stochastic  simulation  in  general,  and  it  is  natural  they 
are  often  jointly  referred  to  by  the  term  “Monte  Carlo.  ”  The  queuing 
problem  is  mentioned  in  many  papers  as  a  natural  example  for  this  method, 
and  is  dealt  with  specifically  in  References  5,  6,  and  7.  Only  a  few  of  the 
ideas  suggested  therein,  however,  can  be  applied  directly  to  the  priority 
queuing  situation.  Some  may  introduce  bias,  others  lose  most  of  their 
efficiency.  The  problem  of  the  distribution  of  the  initial  position  in¬ 
creases  in  dimensionality  and  becomes  almost  too  formidable  to  tackle. 

On  the  other  hand,  there  are  some  advantages  in  these  methods  not  avail¬ 
able  in  the  non-priority  case.  Often,  there  is  at  least  one  estimator 
more  than  there  are  estimands,  which  suggests  the  use  of  a  regression 
model.  Also,  in  some  cases  the  difficulty  is  only  with  the  priority 
discipline.  The  overall  mean  waiting  time  (or  mean  waiting  time  for 
some  of  the  priority  classes)  is  known  analytically,  and  this  knowledge 
can  be  incorporated  in  the  analysis  of  the  experiment  to  reduce  the  var¬ 
iance  of  the  estimator. 

The  utility  of  simulation  depends  not  only  on  getting  estimators 
with  sufficiently  small  variances,  but  also  on  the  ability  to  accurately 
estimate  these  variances.  As  has  been  pointed  out  by  Kahn  and  Marshall, 
“nothing  is  worse  than  thinking  you  have  a  good  estimate  when  in  fact  you 
have  a  very  bad  one."®  The  problem  therefore,  is  twofold:  to  design  the 
experimental  setup  to  yield  practically  no  bias  in  the  estimates  and 
their  confidence  interval,  and  to  reduce  the  variance  of  a  fixed  volume 
of  sampling  without  a  similar  increase  in  computation  labor. 

b.  EXPERIMENTAL  SET  UP  OF  THE  SIMULATION 

The  general  setup  of  the  experiment  followed  closely  the  dis¬ 
cussion  in  Conway.4  The  duration  of  the  run  was  determined  by  the  time 
until  the  departure  of  a  fixed  number  of  arrivals,  i.e.,  “run  until  the 


An  Up-to-date  survey  of  the  development  of  Monte  Carlo  methods  is  presented  in  Reference  5,  where  an 
extensive  bibliography  is  also  to  be  found. 
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first  5000  arrivals  have  departed.'*  Customers  were  grouped  into  sets  of 
fixed  s  i  7.e ,  according  to  the  time  of  their  arrival.  Overall  mean  wait¬ 
ing  time  for  customers  in  a  set  and  the  means  corresponding  to  the  vari¬ 
ous  priority  ('lasses  as  well  as  additional  input  statistics  of  the  set 
were  computed  for  each  of  the  sets,  printed  out,  and  stored  for  the  final 
summary.  Thus,  tin*  individual  sets  are  considered  as  separate  simula¬ 
tions,  the  initial  point  of  winch  is  the  final  position  of  its  predeces¬ 
sor,  and  therefore  the  sets  are  properly  drawn  at  random  from  the 
distribution  corresponding  to  the  simulated  load. 

Experiments  were  started  from  the  “empty  and  idle”  position. 

This  avoids  the  difficulties  involved  in  choosing  from  a  distribution  of 
high  dimensionality,  and  therefore  one  that  would  be  difficult  to  esti¬ 
mate  in  a  preliminary  run.  Beginning  in  this  manner,  however,  introduces 
bias,  and  usually  requires  that  the  first  few  observations  be  chopped  off 
and  not  considered  in  the  summary  analysis.  Theoretically,  the  bias  in¬ 
troduced  by  an  arbitrary  initial  point  (any  point)  can  never  be  totally 
eliminated,  as  indeed  the  “steady  state”  can  never  be  reached  in  a  simu¬ 
lation  run,  no  matter  how  long.  In  practice,  however,  even  mere  chopping 
was  found  unnecessary  for  the  simulation  lengths  used  with  most  of  the 
loading  factors  analyzed.  With  loadings  of  over  0.9  Erlang,  this  prob¬ 
lem  influences  the  estimator  more  strongly.  Statistical  tests  of  the 
transient  effect  with  subsequent  chopping  are  necessary,  as  presented  in 
Appendix  C  of  this  report.  Throughout  the  present  report  the  transient 
effect  is  disregarded. 

It  would  be  proper  at  tins  point  to  define  the  general  setup 
and  basic  output: 

Cast  :  =  nth  arrival 

n 

Setf:=  {Custn:n=  (t-l)m+l . tm }  t  *  1 ,  ...»  T 

p  =  priority  class  p  (all  customers  with  priority  p) 

p  9  1,  ...»  A P  where  p  =  1  is  the  highest  priority 

\P  =  number  of  priority  classes 
Wa i t ( n )  =  Waiting  time  of  Cu  s  t 

Set  =  Set  f|C  =  all  customers  in  Set  that 
t  p  Ml  p  '  t 

have  priority  p  =  1,  ...»  NP 

A  =  number  of  priority  p  customers  in  set,  p  3  1,  ...»  NP 
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w  here 


Cov  [X  t ,  A'  f+  i  ] 
Va  r  X  t 


and  r k 


Cov  \xt,x„k] 
Var 


r  <  1 


with  the  inequality  approaching  equality  for  large  m. 


first  S2 

X 


This  method  involves  at  least  two  substantial  difficulties: 
is  no  longer  an  unbiased  estimator  of  cr 2  but 


ES 2  <  cr2  for  r  >  0 

x  x 

Second,  the  estimate  of  Var  X  depends  on  a  good  estimate  of  r,  and  when 
r  is  large  even  a  small  variance  in  the  estimate  of  r  will  induce  a  large 
variance  in  the  estimate  of  Var  (Af).  Hence,  even  if  independence  is  not 
to  be  assumed,  m  has  to  be  large  enough  to  keep  r  small.  Also,  the  same 
data  may  be  arranged  as  sequences  with  different  values  of  m  (and  T)  and 
the  estimates  compared. 

In  most  cases,  however,  these  elaborate  measures  are  not  justi¬ 
fied,  and  one  may  simply  impose  lower  bounds  on  m  and  7\  with  the  impli¬ 
cation  that  the  volume  of  the  experiment  may  be  determined  not  necessarily 
by  the  variance  of  the  estimates,  but  possibly  by  the  variance  of  the 
estimate  of  these  variances. 

In  the  simulation  of  priority  queues,  we  get  several  estimators 
simultaneously.  The  volume  of  the  experiment  is  determined  by  the 
"worst”  of  these  estimators,  in  terms  of  the  accuracy  required.  Thus, 
although  usually  the  class  with  lowest  priority  will  have  the  highest 
variance,  it  may  be  that  only  a  rough  estimate  is  required  for  that  class. 
The  length  of  the  experiment  would  be  determined  by  a  higher  priority 
class  with  a  smaller  variance.  This  is  especially  true  since  accuracies 
are  usually  specified  as  a  percentage  of  the  estimated  value,  and  for  the 
lowest  priority  not  only  the  variance  but  also  the  estimand  itself  is 
highest. 

Adaptive  rules  to  determine  the  length  of  the  simulation  may 
be  developed,  as  discussed  in  Appendix  B  of  this  report.  The  stopping 
rule  tends  to  be  biased  —  sometimes  as  to  the  estimands  themselves,  but 
more  often  with  respect  to  the  estimate  of  the  variance.  Possibly,  the 
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best  course  is  to  extrapolate  the  required  volume  from  a  short  prelimi¬ 
nary  run.  Note,  however,  whenever  there  is  a  rigid  requirement  on  the 
estimate  of  the  variance  (i.e.,  the  results  have  to  be  presented  as 
P[0.95  It  <  E  [X]  <1.05  »’]  >  0.95),  there  is  a  tendency  to  overestimate 
the  confidence.  This  is  avoided  if  a  certain  variability  in  the  stated 
confidence  is  allowed. 

2 .  Variance  Reduct i on 

a.  INTRODUCTION 

The  basic  principle  of  variance  reduction  methods  involves  tak¬ 
ing  advantage  of  analytical  knowledge  of  part  of  the  simulated  process. 

By  taking  account  of  the  correlation  between  input  and  output,  one  can 
account  for  part  of  the  variance  of  the  output  by  a  proper  analysis  of  the 
input,  or  eliminate  part  of  that  variance  by  a  proper  manipulation  of  the 
input.  To  be  done  effectively,  this  requires  insight  into  the  probabilis¬ 
tic  structure  of  the  simulated  process.  This  is  the  reason  why  “Monte 
Carlo”  often  seems  a  collection  of  “ad-hoc”  methods.  This  is  also  why  the 
application  of  established  principles  to  a  specific  problem  is  not  trivial. 
Even  in  the  limited  field  of  priority  queues,  different  methods  will  have 
different  effectiveness  for  different  loadings. 

The  measure  of  efficiency  of  a  certain  design  relative  to  direct 
simulation  may  be  expressed  as  a  ratio  F,  according  to  Hammersley  and 
llandcomb,  5  where 

“Efficiency  ratio” 

“Labor  rati  on  L  - 

“  Va  r  i  a  n  c  e  ratio”  \p 

where  n  j  and  n  represent  the  computational  volumes  required,  and  where  cr^ 
and  cr2  represent  the  variances  of  the  estimate  for  the  investigated  de¬ 
sign  and  direct  simulation,  respectively.  Whenever  there  is  a  set  of 
estimands,  t p  should  represent  the  ratio  for  the  “worst”  estimate,  in  the 
sense  discussed.  Since  under  different  designs  the  worst  estimate  may  re¬ 
late  to  different  estimands,  it  is  sometimes  preferable  to  consider  the 
e  f  f i c i ency  as 


"  ia\ 


ncr 
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"  1 

V  *  — 

n1 

where  nj  and  n1  are  the  “required”  volumes  of  computation  for  the  speci¬ 
fied  accuracy. 

An  alternative  measure  suggested  by  Ehrenfeld  and  Ben-Tuvia,9 
is  the  relative  reduction  in  variance 


V 


0 


a  “ 


a 


2 
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for  a  fixed  volume  of  computation.  This  is  expressed  as  a  fraction  or 
percentage.  Both  measures  are  used  in  this  report,  and  they  are  either 
directly  specified  or  clear  from  the  context. 

b.  THE  PRIORITY  QUEUE  AS  A  STOCHASTIC  PROCESS 

Considering  a  set  of  size  m  as  an  independent  run,  the  stochas¬ 
tic  process  (  and  subsequently  the  output)  is  defined  by  a  set  of  random 
variables.  The  initial  point  is  a  random  vector,  which  in  the  case  of 
priority  disciplines  has  a  high  dimensionality.  Even  in  the  simplest  case 
of  Poisson  arrivals  and  exponentially  distributed  service  time  the  dimen¬ 
sion  of  the  initial  point  is  AP.  When  m  is  large,  however,  variation  of 
the  output  due  to  the  initial  point  is  fortunately  small,  and  it  is  not 
necessary  to  attempt  to  reduce  this  part  of  the  variance. 

Given  the  initial  point,  the  run  is  determined  by  a  sequence  of 
random  interarrival  times  and  service  times  for  each  of  the  priorities. 
The  values  are  drawn  from  their  respective  distributions.  Alternatively, 
this  can  be  viewed  as  a  single  sequence  of  random  numbers  representing 
interarrival  times,  service  times,  and  the  priority  class.  Variation  in 
waiting  times,  and  subsequently  in  mean  waiting  times,  is  caused  by  vari¬ 
ations  in  the  input  sequence.  These  may  be  classified  as  variations  in 
magnitudes  and  variations  in  the  permutations.  The  magnitudes  are  easier 
to  classify,  analyze  and  control,  but  it  is  the  permutations  that  account 
for  a  larger  part  of  the  variance. 
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CONTROL  VAR  I  ATKS 


r  . 


The  basic  idea  is  as  follows:  if  X  is  an  unbiased  estimate  of 

E  [A]  ,  then 

W  -  A'  -  aV  (where  V  *  V'  -  E  [X] 

V  -  control  variate, 
and  a  is  a  constant.) 

is  also  an  unbiased  estimator.  We  have 


Var  W  -  Var  A  +  a2  Var  V'  -  2a  Cov  [XfV] 


and 


To  find  optimal  value  of  a,  differentiate 


d 

—  Var  H 

9a 


2a  Var  V  -  2  Cov  [XtV]  *  0 


Cov  Lr,V] 
Var  V 


Var  W  -  V  a  r  Af  - 


Cov2  [X,V] 
Va  r  V 


Var  1(1  - 


r2  ) 

*  V 


0  =  1 


(17) 


where  is  the  normalized  correlation  coefficient.  If  Cov  [A\  V']  is 

unknown,  it  can  be  substituted  by  its  estimate  a  W  will  remain  un* 
biased  only  if  cr  is  independent  of  X  and  V .  This  cannot  generally 
be  assumed,  but  the  estimator  is  still  consistent.  With  T  sufficiently 
large  the  bias  can  be  ignored.  Alternatively,  tests  on  bias  may  be  in¬ 
troduced  (see  Appendix  B).  a  can  also  be  estimated  from  a  preliminary 
experiment,  or  a  previous  one  with  a  similar,  though  not  necessarily 
equal,  load. 

The  variance  of  the  estimator  is  estimated  by 


A 
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2 
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A  2 
(Tl 
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V 


(18) 
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When  l7  ‘  and  (cr^  J  are  strongly  and  positionally  correlated, 
this  lias  a  much  smaller  variance  than  t lie  variance  of  a\  In  our  experi¬ 
ments  there  is  a  clear  indication  that  this  was  the  case. 

A  control  variate  is  a  known  random  variable,  t.e.,  EV  and  Var  V 
are  known.  The  idea  can  be  extended  to  a  vector  of  control  variates 
(Tocher  has  presented  this  idea  but  develops  the  CL  by  the  least  square 
me  thod  )  . 10 

Let  V,  E  [V]  ,  C  [V]  be  k  x  1  vectors 

V  =  V  -  E  (V) 


CV]  -  Cov  [X,V] 


Let  2  be  a  k  x  k  covariance  matrix  of  V, 


2 

1  J 


Cov  [V  i%V  j) 


solving  k 
t  i  m  a  1  CL 


r  _2  i  n 
\cv 

Le  t  2 

1 - 

X 

.cv'  \  2  . 

W  =  X  +  aV, 

where  CL 

simultaneous  linear 


,  which  is  also  a  covariance  matrix. 

is  a  k  x  1  vector  of  coefficients.  By 
equations  of  the  gradient,  we  get  for  op- 


a*  *  2*1  •  CV 


and 


Var  If*  -  cr\  -  CVT  2~  'cV  = 

It  follows  from  the  positive  definiteness  of  the  covariance 
matrix11  that 


0  <  Var  w*  <  a2x 

I  will  never  be  singular,  since  it  obviously  would  not  help  us  to  use 
control  variates  which  are  linear  combinations  of  the  others. 


34 


At  tins  point  it  might  help  to  get  more  insight  at  the  struc¬ 
ture  of  multiple  control  variates  by  developing  the  explicit  expression 
for  k  =  2.  I n  this  case, 


r  .  -  r 

X  ,V  1 


v2  v  1  ,  1/  2 


CT 


u  1 


v  \  t  v2 
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^  x  ,  v  2  rxtv\rv\tv2 


v  2 


1  -  r 


t/  1  ,  1/  2 


1  -  TJ 


1  - 


r,,,l  r.,.2  x  .  v  lr  .  ,  v2r  v 1  .  k2 


(19) 


1  -  'V.2 


Suppose  V j  is  the  better  single  control  variate  (r^xl  >  rltV2 
We  want  to  investigate  the  advantage  of  introducing  V 2> 
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We  see  here  that  a  good  “support”  for  a  control  variate  is  not  necessar 
ily  one  that  would  itself  make  a  good  control  variate.  Suppose,  for  ex¬ 
ample,  that  all  three  are  positively  correlated.  By  assumption, 
ri(„i  >  r,tV,  •  r,.„2  maV  be  nearly  as  high  as  r(  but  if  r„1>v2  is 

also  high  there  is  practically  no  advantage  —  e .  g.  , 


r 

1,1/1 

n 

0.8  x 

r 

*  .  v2 
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0.75 

>  ~  ^  (  V\  .  v2  )  ~  V  (  vl  )  “ 

0.  0335 
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exists  when  r  0  and  r  • 
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which  would  be  expected.  Another  and  more  surprising  possibility  is  where 
we  have  a  high  r  .  0 ,  and  r  *  =  0!  Heuristically,  this  means  we  im- 
prove  the  estimate  by  disregarding  that  part  of  the  control  variate  fluc¬ 
tuations  which  have  nothing  to  do  with  it. 

It  is  clear  that  the  more  control  variates  are  operating,  the 
more  difficult  it  would  be  to  find  additional  ones  that  would  still  be 
effective.  But  effectiveness  should  be  considered  on  a  relative  basis, 
in  terms  of  variance  ratio,  and  not  the  efficiency  fraction,  as,  for  ex¬ 
ample,  reducing  the  variance  from  say  10%  to  5%  of  the  original  variance 
involves  reducing  the  necessary  volume  by  half.  Namely: 

1,  .  .  .  ,  v  (fc  )]  =  0.05 

^  [  V  1 . v(*+l)]~  /[l'l,...,v(*)] 

-  -  0.5  (21) 

*  “  ^[  *  1 . V  (  k  +  1  )  ] 

Thus,  control  variates  that  would  initially  be  regarded  as  in¬ 
effective  may  prove  useful  after  the  important  control  variates  have  taken 
effect.  This  is  because  they  either  serve  as  good  complements,  or  because 
they  account  for  a  different  part  of  the  variation.  It  seems  proper  to 
call  them  descriptively  as  “marginal”  control  variates. 

In  searching  for  good  control  variates  in  the  simulation  of 
queuing  systems,  the  sources  of  variation  should  be  kept  in  mind.  The 
problem  is  to  reduce  the  dimensionality  by  finding,  or  composing,  a  few 
control  variates  that  would  maintain  most  of  the  relevant  features  of  the 
high  dimension  input  sequence.  A  part  of  the  relevant  features  is  rep¬ 
resented  by  the  magnitudes  of  the  sequence  elements.  We  expect  waiting 
time  to  be  higher  when  most  service  times  are  high  and  most  interarrival 
times  are  low.  Also,  waiting  times  for  each  priority  would  be  higher, 
the  more  customers  of  the  highest  priority  appear  in  the  sequence.  Over¬ 
all  mean  waiting  time,  however,  in  a  “head  of  the  line”  discipline  would 
remain  unaffected — another  analogy  to  the  professor  who  switched  from 
one  university  to  another  and  raised  the  average  level  in  both  universi¬ 
ties!  A  natural  representation  of  these  attributes  would  be  by  their 
means:  I RT — m  ean  interarrival  time,  SR  V  —  mean  service  time,  and  N  j  — 
observed  number  of  first  priority  customers.  Unfortunately,  but  expec¬ 
tedly,  these  are  all  very  poor  control  variates,  accounting  for  only  a 
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few  percentage  points  of  the  estimators  variance.  The  reason  is  that  the 
variates  do  not  represent  at  all  the  permutations  in  the  sequence,  which 
are  the  main  source  of  variation. 

In  trying  to  account  for  the  permutations,  a  possible  attempt 
would  be  to  match  service  times  and  interarrival  times  in  pairs. 

For  example,  consider  the  recursive  relation  for  single  server 

sy s terns 

W  =  max  (W  .  +  SRVt  .  -  IRT  t  ,  .  ,0) 

j  j  *  *  /  •  i  j  i )  — A 


We  can  make  a  try  by  substituting  W ^  ^  with  a  constant,  and  de¬ 

fining 

V  (Cust  j  )  =  max  {SRV  ,  -  IRV  ]tR) 

j  —  *  )  t  j  ~  * 

Where  R  £  0  is  the  best  estimation  possible,  and  when  P  [W  *  0]  is  sub¬ 
stantial,  R  =  0  is  a  very  convenient  choice.  Now  we  define 

Vt  =  -  V  (Cust  j) 

m  Cust  j  ,  Set { 

and 

Vfp  =  -  £  V  (Cust  j),  where  p  3  1,2  (22) 

<  P 


This  would  not  be  effective  for  the  lower  priority  classes. 

V^,  V2  would  be  the  control  variates  for  the  mean  waiting  times 
of  the  corresponding  priority  classes.  V  can  serve  not  only  for  the  over¬ 
all  mean  waiting  time  but  also  for  all  of  the  priority  classes. 

The  distribution  of  the  V* s  is  easy  to  find;  as  the  mean  of  in¬ 
dependent  and  identically  distributed  random  variables,  the  distribution 
of  each  is  the  truncated  convolution  of  the  distributions  of  service  and 
interarrival  times. 

For  example,  if  arrivals  are  Poisson  and  service  time  is  expo¬ 
nentially  and  identically  distributed,  we  get 
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r  As)  =  1 - '■ —  f '  s  •  o 

1/  T  t  - 

1  +  p 


F  ( s 


1  + 


li  <  s  <  0 


<  A  , 


which  yield 


E  [VI 


1 


pp( 1  +  p) 
1 


[p2  -  (1  -  eA*)] 


Var  [V]  =  - - -  {2(1  ♦  p)[(l  ♦  p3)  -  (1  ♦  \R)eAfl]  -  [p2  -  (1  -  eAfl)]  1  , 

p2p2(l  +  p)2 


and  for  R  =  0 


E  [V\R  =  0] 


P 

P(1  ♦  p) 


Var  [l  fl  =  0] 


P(2  t  p) 
p2(l  ♦  p) 2 


(23) 


The  corresponding  first  moments  for  V  are  immediate  (because 
of  the  independence)  and  for  practical  values  of  mV  t  can  be  considered 
normally  distributed  (this  assumption  is  not  actually  necessary  for  the 
use  of  \  f  as  a  control  variate,  as  only  the  first  moments  are  used).  For 
examp  1 e : 


E  [V'(]  =  E  [V  (Cust)] 

Var  [ V  ]  =  —Var  [ V  (Cust)] 

m 

Generally  speaking,  these  moments  do  not  exist  for  V  ,  as  P[N  =  0]  >  0. 

1  P  tp 

They  do  exist,  however,  for  the  conditional  random  variable  [F  \N  >  0] . 

t  p  t  p 

Clearly, 

E  [V’.PI'V<P  >  °1  =  E  [V  (Cust)] 
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The  derivation  of  l  lie  variance  is  not  computationally  easy.  But  for 
Poisson  arrivals,  the  normal  approximation  with  large  enough  m  can  be 
used  together  with  a  series  expansion  to  get 


Var  lV,p\Ntp  >  0] 


Var  [C(Cust)] 

£  twtp] 


P  r  [Cu  s  t 
£  [N 


Thus  can  be  used  as  a  legitimate  control  variate. 

A  mod  i  f  i  ca  t  i  on  of  K  j  can  be  used  to  obtain  indication  if  priority 
1  customers  arrived  in  groups,  causing  a  substantial  queue  to  be  formed, 
rather  than  “evenly”  distributed.  That  would  be 


Vxl  (Cust  j)  =  max  (SRV)_l  -  IRVjir  0) 


where  j  -  1,  j  are  two  consecutive  arrivals  of  the  priority  class  1  (  re¬ 

gardless  of  any  other  arrivals  in  between). 

In  searching  for  a  good  complement  to  “type  Vn  control  variates, 
note  the  low  correlation  between  waiting  time  and  either  mean  interarrival 
time  or  mean  service  time.  Subsequently,  there  is  practically  no  corre¬ 
lation  whatsoever  with 

SRV  «  eiRT 


where  d  is  a  constant.  Yet,  the  correlation  of  this  compound  variable 
with  V  is  not  eliminated  because  the  influence  of  IRT  in  V  is  truncated 
and  it  can  be  used  as  a  complement  of  the  type  described  above.  There 
is  no  need  to  compute  separately  the  optimal  value  for  8\  rather,  the 
multiple  control  variates  model  can  be  used  with  V,  SRV,  and  IRV ,  and  all 
optimal  coefficients  are  computed  simultaneously.  For  example,  in  a  sim¬ 
ulation  of  a  three  priority,  single  server  exponential  system,  the  use 
of  V  alone  indicated  a  reduction  of  the  variance  to  about  40%  of  its 
initial  direct  simulation  value.  But  combined  with  SRV  and  IRT  the  indi¬ 
cated  reduction  was  down  to  about  8%  of  the  initial  value.  Although  the 
indicated  efficiencies  in  this  case  may  be  s  1  ightly  hi  gher  than  wha  t  could 
be  regularly  expected  with  these  control  variates  (as  would  be  explained 
below)  the  relation  between  the  two  contributions  is  significant. 

A  very  convenient  “marginal”  control  variate  is  /V  f  ^ .  As  noted, 
the  higher  Ar  f j ,  the  higher  would  be  the  expected  mean  waiting  times  for 
each  of  the  priority  classes.  Admittedly  the  correlation  is  low  and  would 
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usually  account  for  only  about  2%  of  the  i  rr  i  t  i  a  l  variance,  but  its  im¬ 
portance  stems  from  the  fact  ,V  (  is  in  many  cases  independent  of  V  (  or 
F  (although  a  somewhat  elusive  dependence  exists  with  V t j j  and  the  ef¬ 
ficiency  is  also  maintained  1  rr  the  margin.  Furthermore,  N  t  j  is  very  con¬ 
venient  computationally  as  a  control  variate  of  X ^ .  Practically  all  of 
the  additional  computation  involved  in  using  a  control  variate  is  in  the 
estimation  of  the  correlation  with  the  estimated  variable,  the  variance 
of  the  control  variate  and  its  covariance  with  other  control  variates  be¬ 
ing  known.  In  the  case  of  ;Y(J  and  X  f,  tins  does  not  involve  even  extra 
multiplications,  as  the  product  X  f  1  N  t  j  is  available  even  before  X  t  ^  [see 
its  Eq.  ( 16)] . 

There  is  another  control  variate  which  is  particular  to  the 
problem  of  priority  queues.  In  many  cases,  the  difficulty  lies  only  with 
the  priority  discipline  (not  with  the  interarrival  and  service  times  sto¬ 
chastic  mechanism).  Analytical  results  are  thus  available  for  the  over¬ 
all  mean  waiting  time  (this  does  not  hold  for  interrupt-repeat  disciplines). 
There  is  a  positive  correlation  between  the  observed  mean  waiting  time 
for  each  of  the  priorities  and  observed  overall  mean  waiting  time.  There¬ 
fore,  A  becomes  a  natural  control  variate. 

It  should  be  pointed  out  here  that  when  X  was  added  to  a  set  of 
control  variates  of  the  types  described  above,  the  increase  in  efficiency 
for  the  high  priority  classes  was  small.  This  means  that  little  is  gained 
by  analyzing  the  overall  process.  It  also  indicated  the  effectiveness  of 
V  as  a  control  variate.  A'  is  most  effective  with  respect  to  a  priority 
class,  or  a  set  of  priority  classes  grouped  together,  that  constitutes  a 
large  percentage  of  the  total  load.  Also,  its  effect  in  the  margin  is 
obvious.  For  the  lowest  priority  class,  A  is  extremely  efficient,  V 

x  *  x  N  P 

often  being  in  the  order  of  0.98. 

A  remarkable  feature  of  control  variates  that  has  been  observed, 
but  not  proved,  is  the  fact  that  not  only  is  the  variance  of  the  estimator 
reduced,  but  also  the  estimate  of  the  variance  is  a  much  better  one.  Since 


A  9 


A  2 
a 

X 


A  j  A 

CV  ‘Z  cv 


this  indicates  that  and  CVl£~lCV  are  positively  correlated  random  var¬ 
iables  and  that  the  correlation  principle  of  the  control  variates  operates 
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also  in  the  estimate  of  the  variance.  As  a  good  estimate  of  the  confi¬ 
dence  interval  is  an  essential  part  of  the  experiment,  this  property  is 
very  significant. 

The  list*  of  the  control  variates  method  may  be  extended  to  var¬ 
iables  of  which  only  the  mean  is  known.  In  this  case  2  has  to  be  esti¬ 
mated  as  well  and  the  problem  of  bias  is  amplified.  Here,  too,  it  can 
be  resolved  by  an  estimation  from  a  preliminary  experiment,  or  extrapo¬ 
lation  from  a  previous  experiment  with  similar  load.  As  long  as  the 
estimate  of  the  optimal  value  for  the  vector  of  coefficients  is  indepen¬ 
dent,  bias  is  not  introduced.  An  error  in  <X  is  thus  never  catastrophic, 
as  at  most,  some  efficiency  will  be  lost.  This  loss  is  of  the  second 
order  only,  because  is  convex  in  (X. 

In  conclusion,  the  control  variates  method  was  found  quite 
effective  for  priority  queues.  The  additional  computation  is  small  and 
just  a  few  carefully  chosen  control  variates  will  eliminate  a  substantial 
part  of  the  initial  variance.  In  a  small  number  of  exploratory  trials 
with  no  more  than  four  control  variates  the  variance  ratio  was  over  2 
under  ” un f a vo r a b 1 e **  circumstances  and  over  10  under  “favorable”  ones. 


d.  STRATIFICATION 

The  basic  idea  is  to  divide  the  observation  space  into  a  set 
of  mutually  exclusive  and  jointly  exhaustive  strata,  depending  on  one  or 
more  “stratification  variables.”  We  must  estimate  that 


by 


-  E  IX\X  e  Sk)  for  each  stratum 


X  k(X  €  S  k)  ,  or  for  any  improved  estimator 


denote 

Pk  =  P[X  e  Sk } 
The  estimator  for  EX  would  be 

It'  =  2  P.WS. 

k 
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V;u  (ill 


a  ml 


1  Pl  Va  r  [IKS.  3 

*  k  * 


Similar  to  what  was  done  before  (for  V'  )  ,  this  is  approximated  by 


Var  [if]  *  2  P  k  Va  r  [X  \  X  €  S  k  ] 


Here  T  (more  important  tlien  m)  has  to  be  large  enough. 

We  see  that  this  is  the  mean  “inter  strata ”  variance,  and  that 
the  “ between  strata”  variance  has  been  eliminated. 

For  a  variable  to  be  eligible  as  a  stratification  variable,  not 
only  the  mean  and  variance  have  to  be  known  in  each  of  the  strata,  but 
also  the  distribution.  With  m  large  enough  this  is  not  a  serious  restric¬ 
tion  because  for  most  of  the  set  statistics  normality  can  be  assumed  (un¬ 
less  autocorrelation  is  exceedingly  high).  The  effectiveness  of  the 
method  depends  on  finding  a  partition  with  the  “between  strata”  variance 
a  large  portion  of  the  total. 

From  the  variables  tested,  none  directly  qualified  in  this  re¬ 
spect.  The  typical  distribution  was  as  follows: 


TA  - 3I»7 -• 

FIG.  IV-1  TYPICAL  DISTRIBUTION  OF  STRATIFICATION  VARIABLE 


We  see  that  most  of  the  variance  actually  lies  in  the  strata 
with  high  V ,  and  that  the  between  strata  variance  is  small.  The  situa 
t l on  is  much  worse  for  “bad”  control  variates  (like  jV j ) . 
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The  method  could  still  be  useful  if  the  number  of  observations 
in  the  important  strata  could  be  artificially  increased  by  manipulating 
the  input.  For  example,  the  simulation  could  be  run  with  higher  load 
than  the  one  investigated,  since  the  weights  P are  analytically  computed 
and  do  not  depend  on  the  number  of  observations.  Unfortunately,  this  can 
not  be  done  under  the  present  experimental  setup  where  the  end  point  of 
one  set  serves  as  an  initial  point  for  the  next.  The  weights  P would 
compensate  for  bias  in  the  input  sequence,  but  not  for  the  bias  in  the 
initial  point. 

Although  some  reduction  in  variance  was  obtained,  on  the  whole 
the  method  was  not  found  too  efficient  for  priority  queues.  Not  only  the 
mean  but  also  the  variance  have  to  be  estimated  in  each  of  the  strata. 

The  reduction  was  usually  not  more  than  just  a  few  percent  (that  could 
not  be  kept  on  the  margin).  It  would  be  more  effective  in  the  simulation 
for  transient  results,  where  the  initial  point  has  to  be  considered  spe¬ 
cifically  anyhow. 

Tli ere  is  one  exception.  In  many  situations  the  probability  of 
having  to  wait  is  known.  We  can  estimate  E  lX\X  >  0]  by,V+,  the  average 
over  positive  values  only,  and  then  estimate  mean  waiting  time  for  each 
priority  class  by 


w  -  p[x  >  oU+ 


which  has  a  smaller  variance  than  X .  For  a  “head-of-the-line"  priority 
queue  discipline  in  an  exponential  system,  the  theoretical  variance  ratio 
would  be  approximately 


0  = 


2 


-  P 


which  was  also  observed  in  actual  experiments. 

This  is  not  a  very  high  reduction,  especially  for  the  more  fre¬ 
quently  encountered  values  of  p.  But  here  again  there  is  prac t i ca 1 ly  no 
extra  work — one  comparison  to  zero  is  added  for  each  customer,  but  the 
number  of  additions  is  reduced.  A  very  slight  increase  is  caused  by  the 
introduction  of  a  counting  index,  altogether  no  change  was  observed  in 
computing  times. 
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The  main  question  is  whether  tins  reduction  would  be  maintained 
(proportionally)  on  the  margin.  In  other  words,  would  the  use  of  control 
variates  be  as  effective  with  respect  to  A  +  as  it  is  for  XI  If  so,  one 
can  obviously  replace  the  estimator  by 

if  =  p[x  >  o]ir+ 

where  If +  is  some  improved  estimator  (though  control  variates)  of  £  [A"  I >  0  ]  . 

It  seems  it  would  be  better  to  use  a  condition  on  the  control  variate  too 
in  this  case,  namely  V +  =  average  of  V  only  over  values  such  that 
V ( Cu s t )  >  R .  The  few  experiments  we  held  indicated  the  control  variates 
method  works  just  as  well  with  X +  as  it  does  with  X  for  the  exponential 
systems,  “head-of-the-lineM  discipline  over  a  wide  range  of  p .  It  may  or 
may  not  hold  for  other  systems  and  other  disciplines. 

A  different  stratification  operator,  also  involving  manipula¬ 
tion  of  the  input,  will  be  discussed  in  connection  with  the  antithetic 
variates. 


d.  ANTITHETIC  VARIATES 

The  basic  idea  is  this.  E[X ]  is  to  be  estimated,  and  Y  is  a 
random  vari  able  such  that  £[V]  =  £[A].  Then 

W  =  1(A  ♦  Y) 

is  an  unbiased  estimate  of  £[A]. 

1  1  1 

Var  If  =  -  Var  X  +  —  Var  Y  +  —  Cov  [A.T] 

4  4  2 

When  Var  A  =  Var  Y 

1  .  .  ,,  1 
Var  W  =  -  {Var  X  *  Cov  [X,  Y) }  =  -  Var  A(1  +  V t y) 

Since  we  have  doubled  the  number  of  observations  (observing  also  Y)  ,  this 
pays  if  V t y  is  negative.  The  is  the  efficiency  fraction. 
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As  pointed  out  in  Maminersley  and  Ilandcomb,5  the  antithetic  var¬ 
iate  method  is  awkward  when  applied  to  problems  of  high  dimensionality, 
as  is  the  case  with  priority  queuing  systems.  One  possible  approach  is 
to  apply  linear  transformations  to  each  of  the  individual  terms  in  the 
input  sequence.  Page  did  this  for  a  no-priority  queuing  system,7  using 
either  the  transformation 


=  1  -  €  T}’  =1-7) 

or  the  transformation 

C  *  *  L 

b  *  V  V  b 

or  both.  denotes  the  random  numbers  that  determine  interarrival  time 
and  rj  the  numbers  that  determine  service  times.  The  efficiency  of  the 
transformation  stems  from  the  fact  that  waiting  time  is  monotone  increas¬ 
ing  with  service  time,  and  monotone  decreases  with  interarrival  time.  In 
turn,  service  time  and  interarrival  time  are  both  monotone  with  the  ran¬ 
dom  numbers  that  generate  them.  The  reported  efficiency  ratios  were  in 
the  order  of  one-half.  In  a  similar  experiment  we  held  for  a  priority 
queuing  system,  the  observed  efficiencies  were  of  the  same  order  for  the 
overall  mean  waiting  time.  The  gain  was  considerably  smaller,  however, 
for  the  individual  priority  classes.  The  reason  may  be  that  the  differ¬ 
ent  runs  are  antithetic  only  with  respect  to  their  **  magn  i  t  ude 99  of  the  in¬ 
put  sequence,  but  not  necessarily  with  respect  to  the  permutations,  which 
account  for  most  of  the  variance.  In  the  priority  case,  the  permutations 
have  an  even  stronger  impact  on  the  results  for  the  individual  classes, 
hence,  the  smaller  gain. 

A  further  transformation  is  presented  in  Hammersley  and 
Handcomb.^  This  consists  basically  of  stratification  of  the  unit  interval. 
The  original  suggestion  was  +  k)/K  k  *  0 ,  1,  ...»  A'-l  ,  but 

this  is  not  directly  applicable  to  the  queuing  problem,  where  in  any 

given  sequence  each  of  the  values  has  to  be  uniformly  distributed  in 
(0,1).  The  modified  transformation  ‘  *  (£  +  k/K )  mod  1,  k  B  0,  1,  ...,  A- 1 

is  legitimate.  This  is  particularly  useful  when  the  function  is  symmetric 

in  the  unit  interval,  which  can  be  aciyeved  by  taking 

f  *  \  [/(£')  +  /ci  -  r )]  • 
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Tin*  natural  way  would  be  to  do  this  separately  for  f  and  T). 

For  each  value  of  the  pair  (k~,k^)  we  may  have  four  runs: 

(1)  (r,V); 

(2)  (l  -r,V); 

(3)  (£\1  -  7J1); 

(4)  (1  -  £,1  -  V> 

This  amounts  to  4A 2  runs. 

Here  a  problem  of  scale  exists.  The  efficiency  of  this  trans¬ 
formation  theoretically  rises  strongly  with  A.  But  we  have  4A"  runs,  the 
length  of  each  is  bounded  from  below.  Therefore,  to  be  very  efficient, 
the  experiment  has  to  be  very  large,  possibly  larger  than  would  be  nec¬ 
essary  by  the  requirements  on  the  results.  For  this  reason,  really  large 
experiments  were  not  attempted. 

An  alternative  approach  rises  from  the  nature  of  pseudo-random 
numbers  used  in  computer  work.  All  computer  random  number  sequences  are 
actually  deterministic,  entirely  dependent  on  the  initial  choice.  Ob¬ 
served  mean  waiting  time  for  each  priority  class  is  therefore  a  positive, 
real,  valued,  function,  albeit  a  complicated  one,  of  a  real  number  in 
(0,1)  —  the  initial  random  number. 

We  can  thus  consider  the  problem  as  a  problem  of  integration 
on  (0,1)  and  avoid  the  high  dimensionality  altogether.  The  antithetic 
transformation  alone  has  no  value,  because  we  now  don’t  have  monotonicity 
but  the  2K  transformations  as  described  above  are  perfectly  adequate.  It 
is  not  obvious  that  the  function  is  continuous,  but  the  congruence  method 
of  random  number  generation  suggests  this  would  be  the  case.  This  alter¬ 
native  approach  was  not  tested  experimentally,  and  might  provide  interest¬ 
ing  basis  for  further  study. 

An  important  question  is  whether  the  use  of  antithetic  (or 
otherwise  dependent)  runs  precludes  the  use  of  other  variance  reduction 
methods,  particularly  control  variates.  One  problem  here  is  that  the 
covariance  matrix  for  the  control  variates  (^)  is  not  easy  to  develop 
analytically  lor  averages  over  dependent  runs.  It  may  sometimes  be  use¬ 
ful  to  choose  as  control  variates  such  functions  of  £  and  rj  that  are  in¬ 
variant  under  the  transformations  used,  even  at  the  cost  of  some  efficiency. 
For  the  four  antithetic  runs  described  above  (without  the  stratification 
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operator),  1/2  -  C  !  and  particularly  1/2  -  rj  |  proved  to  be  reasonable 
approximations.  This  approach  was  not  investigated  further.  Indeed,  one 
of  its  shortcomings  is  that  it  would  involve  substantial  preliminary  ex¬ 
perimentation  to  get  the  proper  insight. 

A  more  direct  approach  is  to  estimate  £  as  well  as  CV  (the  co- 
variance  vector).  Here  the  question  of  bias  becomes  more  acute,  l^ut  jt 
can  again  be  resolved  by  a  preliminary  experiment.  Even  then,  Ct  2~1CV'  is 
a  biased  estimate  of  the  optimal  coefficients,  but  again  it  involves  at 
most  some  reduction  in  efficiency. 

On  the  experiments  carried  out,  the  averages  of  the  control  var¬ 
iates  over  antithetic  runs  did  maintain  the  correlation  with  the  respec¬ 
tive  averages  of  mean  waiting  times.  It  seems,  therefore,  that  both  meth¬ 
ods  can  be  operated  together,  although  not  all  of  the  individual  efficiency 
ratios  would  be  maintained. 

The  antithetic  variates  method  (with  or  without  the  stratifica¬ 
tion  operator)  has  a  certain  appeal  as  an  apparently  general  purpose  method, 
applicable  to  all  kinds  of  systems  and  loads  without  having  to  bother  about 
insight  into  the  stochastic  nature  of  the  process.  This  is  not  always  true. 
Even  if  individual  terms  have  monotone  relationship  with  the  estimated 
variables,  the  structure  for  the  highly  dimensional  input  sequence  as  a 
whole  is  far  more  complicated,  as  indicated  above.  It  may,  under  special 
circumstances,  even  happen  that  the  runs  are  positively  correlated,  which 
would  actually  increase  the  variance. 

e.  SIMULTANEOUS  ESTIMATION  FOB  ALL  CLASSES 

from  a  typical  experiment,  estimates  W  for  each  of  the  priority 
classes  as  well  as  W  for  the  overall  are  available.  To  be  consistent, 
these  estimates  must  satisfy 

jV  P  NP 

W  I  \  =  I  K  W 

p=  1  p  p«  \  p  p 

which  means  that  actually  more  estimates  are  available  than  estimands. 

This  suggests  the  use  of  a  standard  regression  model 
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let  //  be  a  / VP  *  ( /VP  +  1)  matrix 


// 


w  h  e  r  e 


a 


J 


A 

—2. 


NP 

2  \ 

_  ,  i 


J  *  1  . 


/VP 


let  £  be  JVP  x  1  vector  E  =  £  |>  ]  =  £  [X  ] 

P  P  P 

let  If  be  the  ( /VP  +  1)  vector  of  estimators, 


'N, 


If 

P 


JVP+  1 


If 


'■v. 

and  let  W  be  the  covariance  matrix  of  the  estimators 


VW  =  Cov  [if  ,  If  ] 

l  J  l  1  ) 


i,  j  =  1,  2 . NP,  NP  *  1 


Then 


E[W]  =  HE 


Therefore 


if*  =  {HT[VW]-1H}-1HT[VW]-lW 


is  the  best  estimator  of  E  given  If  and 


illT[VW]  -  tf}'1 

is  the  covariance  matrix  for  the  best  estimators. 

The  observed  reduction  in  variance  in  this  step  was  small,  often 
negligible  for  the  higher  priority  classes.  The  additional  work  involves 
estimating  the  elements  Cov  [X  p  x,X  ?  }  and  Cov  [X,Xpi]t  i,  j  =  1,  ...»  /VP, 
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an  operation  the  size  of  which  depends  on  T  and  NP ,  and  inverting  which 
depends  only  on  /VP.  The  method  may  be  more  useful,  therefore,  when  the 
required  accuracy  (and  m )  is  high. 

It  is  often  practical  to  regard  the  estimates  for  the  higher 
classes  as  fixed,  and  employ  the  regression  model  only  for  W N p  and 

NP  -  1  1 

KW  -  2  k  W 

,=  i  1 

This  does  not  reduce  the  work  in  estimation,  but  the  inverted  matrix  is 
2  x  2.  As  the  initial  estimates  for  the  high  priority  classes  are  far 
better,  this  approach  yields  practically  the  same  efficiency  as  employing 
the  la  rge  r  matrix. 

If  the  regression  model  is  discarded,  there  is  a  choice  for  es¬ 
timating  mean  waiting  time  for  the  lowest  priority  with  Yi N p  or 

K  NP-1  k, 

- W  -  1  - —  V 

_  kNP  1 = 1  kNP 

(or  any  convex  combination  of  these  —  the  optimal  one  would  have  been  de¬ 
termined  by  the  regression  model).  The  decision  can  obviously  be  made 
after  the  results  from  the  experiment  are  obtained.  Usually,  when  X  is 
not  one  of  the  control  variates  the  variance  of  W^p  is  much  higher  than 
that  of  the  alternative  estimator.  When  overall  mean  waiting  time  is 
known,  the  variance  of  the  two  estimators  are  of  the  same  order,  but  W ^p 
is  still  the  worse.  Also,  the  estimate  of  Var  l^Np]  should  be  treated 

with  suspicion  in  this  case,  as  small  errors  in  the  estimate  of  r 

r  X'*NP 

result  in  substantial  errors  in  the  estimate  of  Var 

3.  Summary 

The  object  of  this  study  was  to  investigate  efficient  methods  for 
the  simulation  of  queuing  systems  with  priority  disciplines  to  estimate 
mean  waiting  times  under  “steady  state”  conditions.  The  problem  of  esti¬ 
mating  transient  behavior  was  kept  in  the  background.  It  seems,  though,  that 
most  of  the  methods  presented  would  be  applicable  for  the  transient  case 
with  but  little  modification  —  some  of  them  may  indeed  prove  even  more 
efficient. 
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Part  of  the  experiments  were  made  on  the  IBM  7090  computer  at 
Stanford  University,  and  part  on  the  Burroughs  B5500  at  SB I .  For  the 
7090,  the  programming  language  was  SIMSCHIPT,  and  for  the  B5500  an  ALGOL 
simulation  program  already  available  at  SRI  was  used.  Due  to  budget  lim¬ 
itations,  only  a  few  experiments  were  run  for  each  method,  and  the  resu 1 ts 
indicate  orders  of  magnitude  rather  than  precise  conclusions.  As  the  ex¬ 
act  realized  efficiency  depends  on  the  system  to  be  simulated,  the  load 
and  the  specified  accuracy  requirements,  this  seems  the  only  way  for  a 
general  investigation  of  the  nature  of  the  p  re -sent  study. 

The  best  results  were  obtained  with  the  control  variates  method, 
even  when  no  more  than  four  control  variates  were  used  for  each  estimand. 
Conditioning  on  X  >  0  provided  a  further  reduction  in  variance  at  no  ex¬ 
tra  cost.  The  effectiveness  of  the  other  methods  seems  to  depend  more 
on  the  size  of  the  experiment  and  the  required  accuracy.  The  stratifica¬ 
tion  method  provides  the  weights  for  “ importance  sampling"  manipulation 
of  the  input.  That,  too,  would  be  more  applicable  when  required  accuracy 
is  high,  or  for  the  transient  problem,  where  initial  position  lias  to  be 
considered  anyhow. 

The  present  study  was  restricted  to  a  single  service  center  (though 
not  necessarily  a  single  server).  A  natural  extension  would  be  to  queues 
in  a  network.  Obviously,  most  of  the  methods  would  have  to  be  modified, 
and  in  particular  new  control  variates  and  stratification  variables  would 
have  to  be  found.  As  the  stochastic  mechanism  is  more  complicated,  this 
would  probably  be  more  difficult  for  the  network  problem,  and  it  may  be 
expected  that  the  control  variates  method  particularly  would  lose  some  of 
its  relative  advantage  over  the  other  methods.  On  the  other  hand,  the 
idea  of  treating  the  problem  as  a  function  of  the  initial  random  number 
only  gains  some  appeal. 
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\  GUIDES  AM)  PROCEDURES  FOB  THE  APPLICATION 
OF  QUEUING  MODELS 


A  1 NTHODUCT I  ON 

Large  scale  military  information  systems  have  a  number  of  common 
characteristics.  In  these  systems  there  are  usually  one  or  more  inter¬ 
connected  information  processing  centers.  Each  information  processing 
center  is  usually  equipped  with  one  or  more  data  processing  units 
(computers)  and  associated  memory  systems.  The  memory  system  usually 
contains  a  hierarchy  of  memories,  such  as  the  high-speed  core  or  thin- 
film  memory,  the  mass  core  memory,  the  drum  memory,  the  disc  memory  and 
the  tape  memory.  A  typical  hardware  organization  of  the  information 
processing  center  is  as  shown  in  Fig.  V- 1 .  Incoming  messages  arrive  at 


INCOMING 

MESSAGES 


OUTPUT 

DEVICES 


TA-5H7  -  10 


FIG.  V-l  TYPICAL  INFORMATION  PROCESSING  CENTER 

the  input  buffer  where  they  wait  for  the  service  of  the  computer.  The 
incoming  messages  are  usually  categorized  into  classes,  and  the  queuing 
discipline  at  the  input  buffer  is  usually  of  the  priority  type.  An  in¬ 
coming  message  may  require  one  or  more  consecutive  services  from  the 
computer.  The  result  of  each  computer  service  can  be  either  a  new  request 
for  computer  service,  or  a  request  to  retrieve  information  from  the  mass 
memory  or  an  output  message.  An  information  retrieval  request  always 
waits  in  the  mass  memory  buffer.  The  queuing  discipline  at  the  mass 
memory  buffer  can  be  the  priority  type.  The  result  of  an  information  re¬ 
trieval  is  always  a  request  for  the  computer  service.  The  mass  memory  may 
contain  several  nonhomogenous  memory  units.  Access  to  these  memory  units 
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v  a  ii 


be  made  in  parallel.  The  output  messages  queue  up  at  the  output 
buffer  for  transmission  to  the  output  devices. 

The  designer  of  such  an  information  system  is  usually  faced  with 
the  problem  of  determining  the  size  and  configuration  of  the  buffers, 
the  speed  and  the  number  of  computers,  the  service  time  of  the  mass 
memory  and  the  queuing  or  scheduling  discipline  required  when  the  system 
parameters  are  not  known  exactly.  As  a  matter  of  fact,  during  the  early 
stages  of  system  development  the  designer  is  often  interested  in  inves¬ 
tigating  the  expected  performance  and  the  probable  weaknesses  when  the 
parameters  vary  over  a  given  range.  For  example,  during  the  early  stages 
of  system  development  the  computer  softwares  are  not  fully  developed, 
hence  the  computer  service  time  characteristics  are  not  precisely  known. 

It  therefore  is  quite  desirable  to  establish  the  system  response  time  as 
a  function  of  computer  service  time.  The  designer  might  also  be  inter¬ 
ested  in  determining  the  size  of  the  buffers  as  a  function  of  the  size  of 
softwares  and  the  arrival  rate  of  messages.  To  study  the  functional  re¬ 
lationships  between  the  system  parameters  and  the  system  performance  under 
a  specific  system  hardware  and  software  configuration  the  system  designer 
usually  establishes  a  model  of  the  system  and  studies  the  functional  re¬ 
lationships  by  manipulating  the  model. 

The  473-L  Simulation  Model  is  a  particular  case  of  the  generalized 
information  system.  The  procedure  used  in  analyzing  the  473-L  Simulation 
Model  may  therefore  be  used  for  analyzing  the  generalized  information 
system.  This  procedure  is  described  in  the  following  section. 

Although  the  variance  reducing  techniques,  the  control  variates  and 
the  antithetic  variates,  have  not  been  applied  to  the  473-L  Simulation 
Model,  the  possible  application  of  these  techniques  is  indicated. 

B.  GENERAL  PROCEDURES 

The  prerequisite  to  a  successful  analysis  of  any  information  system 
is  the  establishment  of  a  clear  and  precise  description  of  the  system 
from  the  point  of  view  of  servicing  the  incoming  messages.  The  descrip¬ 
tion  should  include  the  following  items  for  each  service  station  within 
a  processing  center: 


52 


(1) 


(2) 


(3) 


(4) 


Input  Message  Characteristics. 


(a) 

The  classification  of  messages. 

(b) 

The  percentage  mix  of  each  class 

of  messages. 

(c) 

The  mean  and  distribution  of  the 
arrival  intervals  of  each  class  o 

inter- 
f  messages. 

(d) 

The  service  requirement  of  each  c 
mess  ages. 

lass  of 

(e) 

The  mean  and  distribution  of  the 
lengths  of  each  class  of  messages 

message 

Buffer  (Queue)  Characteristics. 

(a  ) 

The  organization  of  queues,  i.e.,  is  the  buffer 
divided  into  sections  one  for  each  class  of 

message  ? 

(b) 

The  size  of  queue. 

(c  ) 

The  queuing  (service)  discipline. 

Server  Characteristics. 

(  a  ) 

The  number  and  type  of  servers. 

(b) 

The  mean  and  distribution  of  service  times 

for  each  type  of  service. 

Output  Characteristics. 

(a) 

The  destination  of  the  output  of 

a  server. 

The  service  stations  in  the  typical  processing  center  are  the  com* 
puter  and  the  mass  memory. 


Having  characterized  the  service  centers,  the  operational  measures 
to  be  obtained  from  the  analysis  must  be  defined.  Some  of  the  more  com¬ 
monly  used  operational  measures  are:  the  response  time,  the  queue  length 

and  the  utilization  of  server.  We  usually  like  to  obtain  the  distribution 
of  these  measures,  but  in  practice  this  is  either  too  difficult  or  impos¬ 
sible  to  obtain.  Hence,  we  are  forced  to  settle  for  the  mean  or  an 
extreme  point  of  the  distribution. 


Having  characterized  the  service  centers  and  defined  the  operational 
measures  desired,  the  next  step  is  to  construct  a  model  of  the  processing 
center.  The  analytical  modeling  technique  is  always  preferred  over  the 
numerical  modeling  (simulation)  technique.  Since  the  former  technique  can 
usually  yield  the  analytical  relationship  between  the  system  parameters 
and  the  system  operational  measures,  the  analytical  relationship  can  then 
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be  used  in  system  evaluations  such  as  cost,  effectiveness  studies.  How¬ 
ever,  the  analytical  modeling  technique  is  usually  limited  to  the  modeling 
of  the  simple  service  systems.  Furthermore,  the  results  of  the  analytical 
technique  of  modeling  are  usually  limited  to  the  first  few  moments  of  the 
operational  measures.  Nevertheless,  the  existence  of  even  the  first  mo¬ 
ment,  the  mean  of  the  operational  measure,  may  enable  the  information 
system  designer  to  gain  considerable  understanding  of  the  expected  per¬ 
formance  of  the  system. 

The  available  analytical  results  of  priority -queuing  are  listed  in 
Reference  2.  As  far  as  it  is  known  no  significant  new  results  have  been 
added  since  the  publication  of  this  list.  Any  priority -queuing  service 
center  whose  input,  queue,  and  service  characteristics  do  not  meet  those 
as  listed  cannot  be  analytically  modelled.  However,  there  are  operational 
measures  that  do  not  require  a  knowledge  of  queuing  discipline  even  though 
the  service-center  has  a  priority  queuing  discipline.  For  example,  in  de¬ 
riving  the  mean  utilization  and  the  total  mean  storage  of  a  service  center 
it  is  not  necessary  to  stratify  the  input  into  priority  classes.  Thus  the 
non priority  analytical  results  are  applicable.  If  for  some  reason  the 
storage  is  compartmentalized,  one  compartment  for  each  priority  class, 
then  it  is  necessary  to  derive  the  mean  storage  for  each  priority  class. 

In  such  a  case,  the  priority- queue  results  must  be  used. 

A  typical  information  processing  center  usually  contains  more  than 
one  service  center.  Incoming  messages  usually  require  a  sequence  of 
services  from  these  service  centers.  Each  service  center  usually  contains 
a  buffer  where  messages  may  queue  up.  If  the  buffer  is  small  with  respect 
to  the  expected  loading  of  the  service  center,  blocking  can  occur.  Block¬ 
ing  is  defined  as  the  inability  of  a  service  center  to  service  further 
messages  because  the  buffer  of  the  subsequent  service  center  is  full.  The 
analytical  investigation  of  such  a  service  system  is  quite  formidable  if 
not  impossible.  However,  in  practice  the  buffer  is  usually  quite  large 
and  for  all  intents  and  purposes  the  buffer  can  be  assumed  to  be  infi¬ 
nitely  large.  It  is  thus  possible  to  analyze  each  service  center 
independently. 

frequently  the  information  system  is  too  complex  for  analytical 
modeling.  Thus,  we  must  resort  to  the  simulation  modeling  technique.  As 
was  pointed  by  Conwa  y  ,4  t  lie  re  are  three  phases  in  an  investigation  by  simu¬ 
lation  that  take  place  after  the  problem  has  been  identified  and  a  model 
f  ormul a  ted : 
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(1)  Model  implementation — description  in  a  language  acceptable 
To  tTTe  a  p  p  r  o  p  r  i  a  t  e  computer  . 

(2)  Strategic  planning — design  of  experiment  that  will  yield 
t  lie  desired  inform  a  t  i  o  n  . 

(3)  Tactical  planning — determination  of  how  each  of  the  test 
runs  specified  in  the  experimental  design  is  to  be 

exec  u  t  ed . 

There  exists  a  number  of  s i mu  1  a t i on - o r i en ted  computing  languages, 
among  them  are  the  SIMSCRIPT,  the  General  Purpose  System  Simulation 
(GPSS1 I ) ,  the  Control  and  Simulation  Language  (CSL)  and  the  MILITRAN. 

No  attempt  has  been  made  to  evaluate  the  suitability  of  these  languages 
to  the  simulation  of  information  systems.  In  this  project  the  original 
simulation  model  was  programmed  in  ALGOL,  since  it  was  felt  the  model 
under  study  was  relatively  simple  and  would  be  more  efficient  to  use 
ALGOL.  During  the  strategic  planning  and  tactical  planning  phases,  how¬ 
ever,  it  was  necessary  to  make  frequent  changes  in  the  sampling  schemes. 

An  ALGOL  simulation  program  is  not  the  easiest  program  for  extensive  modi¬ 
fications.  Thus,  the  simulation  model  was  reprogrammed  in  SIMSCRIPT  during 
the  later  part  of  this  project.  The  reason  for  choosing  SIMSCRIPT  over  the 
other  simulation  languages  are: 

(1)  SLMSCRI PI  is  the  most  commonly  used  simulation  language,  and 

(2)  a  SIMSCRIPT  compiler  is  available  with  the  Stanford's  7090 
c  omputer . 

An  aspect  of  the  experimental  design  for  simulation  experiment  is 
reduction  of  the  variance  of  a  fixed  volume  of  sampling  without  corre¬ 
sponding  increase  in  computation  labor.  Three  variance  reducing  techniques 
have  been  discussed  in  Sec.  IV-D,  namely,  stratification,  control  variates 
and  antithetic  variates. 

The  stratification  technique  requires  a  knowledge  of  the  distribution 
of  the  stratification  variable.  Usually  the  distribution  is  not  known. 
Although  it  may  be  possible  to  assume  a  normal  distribution,  it  is  not 
possible  to  estimate  the  confidence  attributable  to  the  results.  This 
technique  is  thus  deemed  impractical  for  information  system  simulation. 

There  does  not  seem  to  be  any  set  procedure  in  the  application  of  the 
control  variates  technique  to  the  simulation  of  information  systems.  It 
usually  involves  a  careful  examination  of  the  probabilities  structure  of 
the  simulated  process.  The  purpose  of  the  examination  is  to  determine  the 
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spec i lie  input  or  the  combination  ol  inputs — say  V — in  the  simulation 
that  clearly  contributes  to  the  variation  in  the  output.  Having  found 
V',  and  knowing  the  distribution  of  the  input(s)  that  make  up  the  V ,  it 
is  possible  to  derive  the  E[V ].  Thus,  the  estimate  of  the  output  is 
H  =  A  -  ft  1  V  -  E[V])  where 


X  is  the  simulated  estimate  of  the  output 

V  is  the  simulated  estimate  of  the  control 
variate  V 

E  1j  is  the  theoretical  mean  of  the  control 
variate  V 


Ot  1  s  a  c  onst  ant . 


A  good  control  variate  for  a  single  exponential  server  priority 
queuing  system  has  been  found  to  be  (Sec.  IV-D): 


where 


t  p 


1 

P  Cuat  I  €  Set 

1  t  p 


k( Cus  t  j ) 


V(Cust  j ) 
SEE  . 

j  -  1 


max  (SEBj  _  |  -  IRV} 

tlie  service  time  of  the  j  -  1-th  customer 


IRB 

j 


-  i 


the  time  interval  between  the  (j  -  1 ) - 1 h 
arrival  and  the  j t h  arrival 


R  -  constant 


Other  simple  control  variates  are:  interarrival  time,  service  time, 

number  of  highest  priority  customer  and  the  overall  mean  waiting  time. 

\Se  have  found  it  possible  to  use  more  than  one  control  variate. 


In  a  more  complex  service  system,  say  multiple  parallel  servers,  a 
good  control  variate  may  be  rather  difficult  to  find.  Thus,  we  may  have 
to  use  the  simple  control  variates  and  be  contented  with  lower  reduction 
in  variance.  for  example,  in  the  case  of  the  473-L  Simulation  Model  the 
distribution  of 
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(1)  t  lie  i  ii  i  <*  rn  r  r  i  vu  I  interval  of  requests, 

(2)  the  service  time  at  the  computer, 

(3)  the  service  time  at  the  disc,  and 

(4)  the  number  of  disc  searches, 

can  be  used  as  the  control  variates  for  the  mean  response  time  measurement. 

While  as  the  distributions  of 

(1)  the  message  types, 

(2)  the  request  message  length, 

(3)  the  interarrival  intervals  of  requests, 

(4)  the  search  type, 

(5)  the  data  length, 

(6)  the  program  length,  and 

(7)  the  display  message  length, 

can  be  used  as  the  control  variates  for  the  mean  storage  size  measurement. 

Our  investigation  of  antithetic  variates  technique  has  not  advanced 
far  enough  to  allow  one  to  suggest  any  specific  procedure  for  the  appli¬ 
cation  of  this  technique  to  information  system  simulation.  However,  the 
procedure  outlined  in  Tocher9  may  be  used,  if  there  are  k  distributions 
flp  ft2>  •••,  ft  k  to  be  used  as  sources  of  antithetic  variates.  In  the 
473- L  Simulation  Model,  k  -  4  was  used  for  the  mean  response  time  measure¬ 
ment  and  k  -  7  for  the  mean  storage  size  measurement.  It  has  been  recom¬ 
mended  that  k  should  be  restricted  to  2  or  3.  Thus,  we  must  select  those 
distributions  that  affect  the  measurement  most  directly.  In  the  case  of 
mean  storage  size  measurement,  the  interarrival  intervals  of  requests, 
the  search  type,  and  the  program  length  distributions  are  the  three  dis¬ 
tributions  that  are  most  directly  correlated  to  the  measurement.  There¬ 
fore  these  three  distributions  should  be  used  as  the  primary  control 
variables.  Choose  a  sample  size  n  m  2 * .  In  each  of  n  sets  of  2*  factorial 
experiment.  Associate  a  vector  V^,  V 2#  ...,  V with  each  sample.  If 

V  -  0,  use  the  variable  <f  in  forming  /?  ;  rf  V  -  1  use  the  variable 
1  -  £ .  in  the  transformation  to  ft  . 

i  i 

For  k  -  3,  the  table  which  follows  illustrates  the  arrangement. 
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Table  VI 


A  DKSIGN  OF  AN  KXPKHIMKNT 
IOH  ANTI  TIIKTIC  VAHIATKS  SIMULATION 


DISTRIUI1TI0N 

variable 

SAMPLE 

1 

2 

3 

4 

5 

6 

7 

8 

». 

1 

0 

1 

0 

1 

0 

1 

0 

1 

«  2 

2 

0 

0 

1 

1 

0 

0 

1 

1 
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Associated  with  /?  j  there  are  four  antithetic  pairs,  namely,  (1,2),  (3,4), 
(5,b),  and  (7,8).  Let  their  averages  be  denoted  1 ' ,  2 ' ,  3 ' ,  and  4 ' . 

The  pairs  (\* ,2* )  and  (3,,4#)  are  antithetic  pairs  from  distribution 
ft  2  and  their  averages  form  an  antithetic  pair  for  ftv  Let  z^k  be  the 
sample  value  of  the  jth  sample  of  the  k th  set,  where  j  =  1,  2,  ...»  8  and 
k  -  1,  2,  .  Let  Z  =  (1/8)  2  z  .  .  Then  the  estimate  is  _L 

_  *  j  J  k  —  .  2 

Z  -  (  J  /n )  1  Z,  .  Tlie  estimated  sample  error  is  found  from  [2  (Z,  -  Z)2] 

k  k  * 

and  the  expected  variability  of  the  response  is  derived  from 

j_ 

(1/8)  I  [I  (  z  -  7.)2]  2 

j  »  J 

where 


k 


For  those  who  are  interested  in  the  design  of  simulation  experiment 
and  the  tactical  planning  of  simulation  experiment,  References  3  and  9 
should  be  consulted.  Tocher,9  states  the  whole  field  of  experimental 
design  for  simulation  experiments  is  in  its  infancy  and  offers  a  fertile 
field  for  further  research.  Therefore,  more  efficient  procedures  will 
certainly  be  developed  in  both  the  application  of  variance  reducing  tech¬ 
niques  and  in  the  tactical  planning  of  simulation  experiments. 
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VI  CONCLUSIONS  AND  IIECOMMENDATIONS 


A  CONCLUSIONS 

Often  it  is  possible  to  idealize  the  relatively  complex  network  of 
service  centers  of  an  information  system  into  independent  service  centers. 
B\  assuming  these  to  be  exponential  service  centers  it  is  possible  to 
apply  the  existing  queuing  theory  to  model  the  information  system.  Al¬ 
though  the  operational  measures  obtained  from  the  analytical  modeling 
may  not  be  exact,  nevertheless  they  provide  an  order  of  magnitude  esti¬ 
mate  of  the  system  performance,  and  also  a  check  on  the  Simulation  Model 
is  obtained. 

Analytical  modeling  techniques  are  usually  limited  to  the  analysis 
of  relatively  simple  queuing  systems.  Even  then  the  analytical  technique 
can,  in  most  cases,  yield  only  the  mean  of  the  measure.  The  distribution 
of  the  measure  is  usually  left  in  the  form  of  the  transformation  of  the 
generating  function.  It  would  seem  that  a  fruitful  research  effort  is 
to  deduce  the  approximate  distribution  from  the  transform  of  the  genera¬ 
ting  function  without  actually  carrying  out  the  transform. 

For  more  complex  queuing  systems  we  must  rely  upon  the  simulation 
techniques  for  the  analysis  of  the  system  performances.  The  principal 
drawback  of  the  simulation  techniques  are 

(1)  the  amount  of  effort  involved  in  the  construction  and 
check  out  of  the  simulation  model,  and 

(2)  the  relatively  large  sample  size  required  to  obtain 
the  operational  measures  that  are  sufficiently 
accurate. 

With  the  appearance  of  the  s i mu  1  a t i on - o r i en t ed  computer  languages,  part 
of  the  first  drawback  has  been  overcome.  The  sample  size  required  in 
the  simulation  can  be  drastically  reduced  by  careful  planning  of  the 
simulation  strategies  and  tactics.  The  whole  field  of  simulation  strat¬ 
egies  and  tactics  is  still  in  its  infancy  and  offers  a  fertile  field  for 
further  research. 
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B.  RKCOMMKNDAT I ONS 


It  seems  that  the  simulation  modeling  technique  is  going  to  be  the 
principal  technique  for  analyzing  the  queuing  situations  in  information 
systems,  unless  there  is  a  major  breakthrough  in  the  analytical  technique 
which  is  not  very  likely.  Therefore,  it  is  recommended  that  future 
research  effort  in  the  application  of  queuing  theory  to  information 
systems  design  be  directed  to  the  improvement  of  the  simulation  efficiency. 
In  particular,  the  application  of  the  control  variates  and  the  antithetic 
variates  method  in  the  simulation  of  a  network  of  service  centers  should 
be  investigated. 


The  derivation  of 
Lion  of  the  generating 
although  it  is  not  yet 


the  approximate  distribution  from  the  transforma- 
f unction  should  be  a  worthwhile  research  effort, 
known  how  this  can  be  accomplished. 
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APPENDIX  A 


QUEUING  ANALYSIS  OF  THE  SYSTEM 


I  INTRODUCTION 

The  principal  purpose  of  the  queuing  analysis  of  the  47  3-L  system 
is  to  investigate  the  possibility  of  applying  the  queuing  theory  to  ob¬ 
tain  some  of  the  queuing  measures  of  a  real  information  processing  system. 
The  473-L  system  model  has  been  chosen  for  this  purpose  because  a  simula¬ 
tion  model  of  this  system  has  been  constructed  and  simulated  results  of 
this  system  exist.  Thus,  it  is  possible  to  check  the  analytical  results 
with  the  simulated  results. 

The  473-L  system  model  is  extracted  from  the  monthly  status  report 
previously  cited.1  In  this  analysis,  the  analytical  expressions  for  the 
utility  of  the  computer,  the  utility  of  the  disc,  the  mean  response  time, 
and  the  mean  storage  size  were  obtained. 

II  DESCRIPTION  OF  473-L  SIMULATION  MODEL 

Figure  A- 1  shows  the  473-L  Simulation  Model.  Message  requests  from 
the  “Source  of  Data”  arrive  at  Buffer  1,  Q- 1 ,  in  random  time — that  is, 
the  interarrival  time  is  exponentially  distributed  with  mean  arrival  rate 
of  For  simplicity,  message  requests  shall  be  called  new  requests,  /?. 

The  new  requests  are  categorized  according  to  their  priority  level  and 
message  type.  There  are  three  priority  levels  and  three  message  types; 
thus  R  .  is  the  ith  priority  level  and  j  th  message  type  request,  where 
i  *  1,  2,  3  and  j  *  1,  2,  3.  In  Reference  1,  the  j  =  1  message  type  is 
called  the  electronic  display  type  (ET),  the  j  =  2  message  type  is  called 
the  multicolor  display  type  (MC),  and  the  j  -  3  message  type  is  called 
the  console  printout  type  (CP).  The  i  *  1  request  is  the  highest  pri¬ 
ority  request  and  i  =  3  is  the  lowest  priority  request.  All  arriving 
requests  are  stored  in  Q- 1 ,  which  is  a  part  of  computer  core. 

New  requests  will  always  interrupt  the  computer  operation.  If  the 
new  request  is  of  lower  or  equal  priority  to  the  one  that  the  computer 
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(1)  SEARCH  TYPE 
$0%  OAT  A 

GAUSSIAN  MEAN  *  16  WOROS,  SIGMA  *  3  WORDS 
50%  PROGRAM. 

Gaussian-  mean*,  sigma* 

(2)  DISPLAY  LENGTH: 

ET  UNIFORM  1  -297  WOROS 

MC  GAUSSIAN  MEAN  *  2300  WORDS ,  SIGMA  »  700  WORDS 
CP  GAUSSIAN  MEAN  *297  WORDS  SlGMA  *  90  WORDS 

(3)  SERVICE  TIME 

FOR  SEARCHES  EXPONENTIAL  MEAN* 

FOR  STORAGE  LENGTH  IN  WORDS  X  072 


(1)  COMPUTER  CORE 

(2)  UNLIMITED  SIZE 

(3)  FIFO  DISCIPLINE 


TD  POP  UP  LIST 


SOURCE 

(1)  TIME  BETWEEN  ARRIVALS 

EXPONENTIAL  MEAN* 

(2)  THREE  LEVEL  PRIORITY 

63  2%  l  OW 
26.3%  MEOIUM 
10.5%  HIGH 

i3)  display  or  message  type 

80%  ET  OR  L 119 
10%  MC 
10%  CP 

(4)  REOUEST  LENGTH 
50%  2  WORDS 

50%  UNIFORM  1  -297  WORDS 
ft  RESPONSE  TIME  IS  TIME  DELAY  FROM  POINT  A  TD  POINT  B 

parameters  varied 

(1)  MEAN  TIME  BETWEEN  ARRIVALS:  7.  3,  4,  5,  6,  7,  9,  12  SEC 

(2)  COMPUTER  SERVICE  TIME-  40,  50.  65,  75  80,  90,  180,  200  220  mec 

(3)  SEARCHES:  MEAN  =  9,  1$,  SlGMA*  3,  S 

(4)  PROGRAM  LENGTH  .MEAN  *  1000,  2500  WORDS.  S IGMA  *  300,  600 

(5)  DISC  SERVICE  TIME  90,  180,  360,  540 


BUFFER. 

(1)  COMPUTER  CORE 

(2)  UNLIMITED  SIZE 

(3)  POP -UP  LIST  LENGTH 

(o)  FDR  TASK  IN  PROCESS 

*  SEARCHES  X  S  WORDS 
(b)  FOR  REOUESTS  INOUEUE 

#  REOUESTS  X  1  WORD 

(4)  PRIORITY  ORDERING  IN  BUFFER 


COMPUTER 


am 


(1)  SERVICE  TIME 

EXPONENTIAL:  MEAN* 

(2)  SEARCHES 

25%  0 

7S%  GAUSSIAN:  MEAN* 
SIGMA* 

(3)  ALL  BUFFER  INPUTS  CAUSE 
HOUSEKEEPING  INTERRUPT 
TIME  1  ms*c 


READOUT  RATE 
ET  *  0 

MC*  S  6  ms*c  X  I  WORDS 
CP  =  560  m»ec  X  *  WOROS 


T  •  - 


mr*i 


FIG.  A-l  473-L  SIMULATION  MODEL 
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is  working  on  when  interrupted,  the  new  request  is  routed  to  Buffer  2, 

Q-2,  which  is  the  other  part  of  compute!  core.  This  new  request  is 
labelled  as  a  pop-up  request  and  will  wait  in  Q-2  until  it  is  read  into 
the  disc,  under  the  FIFO  (First  in,  First  out)  discipline.  At  the  same 
time,  a  pop-up  entry  is  added  to  the  pop-up  list  in  Q- 1 .  If  the  new 

request  is  of  higher  priority  to  the  interrupted  message,  the  computer 

will  analyze  the  request  and  initiate  the  first  disc  search,  if  any, 
required  to  fulfill  the  request.  In  all  cases  of  interruption,  the  in¬ 
terrupted  request  is  put  hack  in  Q-l,  at  the  head  of  its  priority  queue. 

Twenty-five  percent  of  the  requests  require  no  disc  search,  they  are 
the  FT  type  requests.  These  requests  are  denoted  by  The  remaining 

requests  have  a  disc  search  distribution  of  the  Gaussian  type  with  mean  n 
and  variance  crrf.  These  requests  are  denoted  by  R  t  j  n-  The  R .  . Q  has  higher 

priority  than  Rljn  f°r  1  -  1,  2,  3.  The  R ^  . n  requires  one  additional 

search  to  retrieve  intermediate  answers  and  format  the  reply  message. 

The  new  request  whose  first  search  has  been  initiated  is  placed  at 
the  end  of  its  priority  queue  in  Q-l,  and  this  request  is  labelled  as  the 
"repeat”  request.  The  computer  will  accept  the  first  request  in  the 
highest  priority  queue  in  Q-l  as  the  next  request  to  be  analyzed.  When 
the  computer  accepts  a  14  repeat”  request  that  requires  further  disc  search, 
it  will  initiate  the  next  disc  search  and  place  the  request  at  the  end  of 
its  priority  queue  in  Q-l.  The  label  of  this  request  remains  as  the 
“repeat”  request.  If  no  further  disc  search  is  required,  the  computer 
will  initiate  the  last  search  and  label  this  request  as  the  “computer” 
request  and  place  the  request  at  the  end  of  its  priority  queue  in  Q-l. 

When  the  computer  accepts  a  pop-up  entry,  the  computer  initiates  a  disc 
search  to  retrieve  the  pop-up  request.  Thus,  a  pop-up  request  requires 
one  additional  disc  search. 

Ihe  disc  search  requests  and  the  pop-up  store  requests  are  queued  up 
at  Q-2  in  the  order  of  arrival.  When  a  search  is  completed,  the  request, 
a  44  repeat  ”  request,  is  routed  to  Q-l.  If  the  "repeat”  request  has  higher 
priority  than  that  of  the  request  being  processed  in  the  computer,  the 
computer  process  is  interrupted  and  the  request  in  the  computer  is  removed 
from  the  computer  and  is  placed  in  Q-l.  The  processing  of  the  "repeat” 
request  is  started.  Otherwise,  the  “ repeat  ”  request  is  placed  at  the  end 
of  its  priority  queue  in  Q-l. 
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The  " repeat”  request  returning  from  the  disc  is  first  analyzed  by 
the  computer.  If  it  requires  further  search,  then  a  disc  search  request 
is  generated  and  placed  in  Q- 2.  If  it  requires  no  further  search,  this 
means  a  reply  message  is  ready  for  transmission;  the  reply  message  is 
transmitted  to  the  output  devices,  representing  both  displays  and  Lll9s. 

The  time  when  the  message  request  arrived  at  Q- 1  to  the  time  when  a 
reply  message  is  ready  for  transmission  is  defined  as  the  response  time. 

If  the  reply  message  is  the  MC  or  CP  type,  it  has  additional  transmission 
delay  that  is  added  to  response  time  to  give  turn-around  time. 

Ill  SYSTEM  PARAMETERS 

A.  Input 

As  mentioned  before,  the  interarrival  time  of  new  requests  is  as¬ 
sumed  to  be  exponentially  distributed  with  mean  arrival  rate  of  A.  The 
new  requests  are  categorized  according  to  their  priority  level  and  their 
message  type.  The  percentage  of  requests  that  are  the  ith  priority  level 
is  ,  where  Y j  *  10.5%,  Y2  =  26.3%,  and  Y3  =  63.2%.  The  percentage  of 
requests  that  are  the  j  th  type  is  ,  where  H'j  =  80%,  Yi 2  =  10%,  and 
V  e  10%.  Fifty  percent  of  the  new  requests  have  a  length  of  two  words; 
the  other  50%  have  a  length  distributed  uniformly  between  1  and  297  words. 

B.  Buffer  (Q-l) 

Q- 1  is  part  of  the  computer  core,  the  size  of  which  we  assume  to  be 
unlimited.  Q- 1  is  used  to  store:  (a)  the  new  requests  from  the  source 

of  data,  (b)  the  pop-up  list,  and  (c)  the  results  of  disc  search. 

The  pop-up  list  contains  two  types  of  request.  The  first  type  of 
request  are  those  requests  whose  service  has  been  initiated  by  the  com¬ 
puter;  the  length  of  this  type  of  message  is  the  number  of  searches  times 
five  words.  The  second  type  of  request  are  those  requests  that  have 
been  sent  to  the  disc  for  temporary  storage.  Their  length  is  one  word. 

The  second  type  of  request  becomes  the  first  type  of  request  when  it  has 
been  retrieved  from  the  disc  for  processing. 

The  length  of  the  request  returning  from  the  disc  search  depends 
upon  the  type  of  disc  search.  If  the  search  is  for  retrieving  the  re¬ 
quests  in  the  pop-up  list,  the  length  distribution  is  the  same  as  the  new 
requests.  If  the  search  is  for  retrieving  data,  50%  of  the  searches  are 
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of  this  type,  the  length  is  Gaussian  distributed  with  mean  Mds  *  16  words 
and  standard  deviation  o,  =3  words.  The  other  50%  of  the  searches  are 

a  s 

the  program  retrieval  type;  their  length  is  also  Gaussian  distributed 
with  mean  M  and  standard  deviation  cr  Finally,  the  last  search,  the 

p  s  pi 

search  to  retrieve  the  intermediate  answers  and  format  the  reply  message, 
has  a  Length  distribution  that  depends  upon  the  type  of  request:  for  the 

FT  type  the  length  is  uniformly  distributed  from  1  word  to  297  words; 
for  the  MC  and  CP  types,  the  length  has  a  Gaussian  distribution  of  mean 
l2  ~  2300  words  and  Z3  *  297  words  and  standard  deviation  cr  2  *  700  words 
and  cr  3  =  90  words,  respectively. 

All  requests  in  Q- 1  are  arranged  according  to  the  priority  level  and 
in  the  order  of  the  time  of  arrival  at  Q- 1  within  each  priority  level. 

C.  Computer 

The  service  time  of  the  computer — the  length  of  time  required  to 
analyze  a  request  and  initiate  a  disc  search  —  is  exponentially  distrib¬ 
uted  with  mean  h c .  Twenty-five  percent  of  the  requests  require  no 
search;  the  other  75%  have  a  Gaussian  distribution  of  disc  search  with 
mean  n  and  standard  deviation  cr  . 

s 

Fach  interruption  of  the  computer  operation  consumes  one  msec  of  the 
computer  time  for  the  h o use keep i ng  function.  The  housekeeping  time  is 
symbolically  represented  by 

U.  BuFFEn  (Q-2) 

The  inputs  to  Q- 2  are  part  of  the  output  of  the  computer.  The  input 
consists  of  three  types  of  request:  the  first  is  the  storage  type  R  , 

the  second  is  the  pop-up  request  retrieval  type  Br ,  and  the  third  is  the 
disc  information  search  type  Bd.  The  message  length  of  these  three  types 
of  requests  are  as  follows:  B  has  the  same  distribution  as  the  new  re¬ 
quest  to  Q- 1 ,  B r  is  one  word  long,  and  B d  is  5  words  long. 

The  buffer  size  is  assumed  to  be  unlimited.  Bequests  are  arranged 
according  to  their  order  of  arrival  at  Q-  2. 

E.  Disc 

The  disc  is  assumed  to  have  unlimited  size.  The  service  time  for 
storing  a  pop-up  request  hdg  is  proportional  to  the  length  of  the  request 
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or  0.072  msec  x  length.  Hie  service  time  for  conducting  a  disc  search  is 
exponentially  distributed  with  a  mean  hd. 

The  output  of  the  disc  search  is  part  of  the  input  to  Q-  1. 

IV  ANALYTICAL  ANALYSIS 

The  473-L  Model  is  basically  a  tandem  queue  service  system.  The 
customers,  new  requests,  may  require  successive  services  at  the  two 
servers,  the  computer  and  the  disc.  for  example,  a  new  request  that  re¬ 
quires  n  disc  searches  must  obtain  service  first  at  the  computer  and  then 
at  the  disc,  and  then  repeat  this  service  sequence  for  n  +  1  times  before 
its  service  requirement  is  completed.  Since  the  size  of  Q- 1  and  Q- 2  has 
been  assumed  to  be  unlimited,  the  Q-  1  and  Q- 2  are  said  to  be  statisti¬ 
cally  independent.  A  request  with  n  disc  searches  is  said  to  have  passed 
through  2(n  +  1)  s t a t  i  s t i c a  1 1 y - i ndependen t  stages  of  queue. 

Throughout  this  analysis,  we  assume  that  the  requests  arrive  at  the 
Q- 1  and  Q- 2  according  to  Poisson  distribution.  This  assumption  is  justi¬ 
fied  by  the  following  argument.  The  new  requests  arrive  at  the  Q- 1  ac¬ 
cording  to  Poisson  distribution.  The  service  time  distribution  of  the 
computer  is  exponential.  It  has  been  shown  in  Burke,^  that  the  inter- 
departure  of  a  Poisson  input  and  exponential  service  system  is  exponential. 
Thus,  the  interdeparture  of  requests  from  the  computer  is  exponential. 
Since  the  inputs  to  the  Q- 2  are  the  outputs  of  the  computer,  the  input  at 
the  Q- 2  is  Poisson.  Now,  the  service  time  of  the  disc  is  assumed  to  be 
exponential;*  thus,  the  interdeparture  of  requests  from  the  disc  is  expo¬ 
nential.  It  follows  that  the  return  requests  from  the  disc  will  arrive 
at  the  Q- 1  according  to  Poisson  distribution. 


The  service  time  distribution  of  the  disc  is  strictly  speaking  not  Poisson.  No  more  than  7%  of  the 
requests,  depending  upon  the  utility  of  the  disc,  are  of  the  storage  type.  Fifty  percent  of  this  type 
of  request  requires  s  constant  service  time  of  0.H4  msec,  and  the  other  50%  hss  s  uniform  service 
time  distribution  from  0.072  msec  to  21.38  msec.  The  other  93  or  more  percent  of  the  requests  hss  an 
exponential  service  time  distribution  with  mean  no  less  than  90  msec.  Since  the  percentage  of  store 
request  is  relatively  small,  its  influence  on  the  overall  service  time  distribution  is  negligible. 
Hence,  the  exponential  assumption  is  s  good  first  approximation. 
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The  following  notations  are  used  in  the  analysis: 


\  c 


i  ;  n 


Y  - 


W  = 
) 


M 


dt 


a 


M 


p  * 


Or 

P  * 


1 

) 


a 

) 


h 


X  * 


n  - 


a  - 

s 


/? 

I 

/? 

r 


/l  * 

s 


The  mean  arrival  rate  of  new  requests. 

The  ith  priority  and  j th  type  new  request  with 
n  disc  searches,  i  =  1,  2,  3  and  j  *  1,  2,  3. 

The  percentage  of  new  requests  that  are 
ith  priority  requests. 

The  percentage  of  new  requests  that  are  j th  type 
requests. 

The  mean  of  the  Gaussian  distribution  of  the 
retrieved  length  of  the  data  type  disc  search. 

The  standard  deviation  of  the  Gaussian  distribu¬ 
tion  of  the  retrieved  length  of  the  data  type 
disc  sea  rch . 

The  mean  of  the  Gaussian  distribution  of  the  re¬ 
trieved  length  of  the  program  type  disc  search. 

The  standard  deviation  of  the  Gaussian  distribu¬ 
tion  of  the  retrieved  length  of  the  program  type 
disc  search. 

The  mean  length  of  the  last  disc  search  for  the 
j t h  type  requests. 

The  standard  deviation  of  the  length  of  the  last 
disc  search  for  the  jth  type  requests. 

The  mean  service  time  of  the  computer. 

The  percentage  of  new  requests  that  require  no 
disc  search. 

The  mean  of  the  Gaussian  distribution  of  the 
number  of  disc  searches. 

The  standard  deviation  of  the  Gaussian  distribu¬ 
tion  of  the  number  of  disc  searches. 

The  interrupt  housekeeping  time. 

The  storage  type  requests  at  Q-  2. 

The  pop-up  request  retrieval  type  at  Q-2. 

The  disc  data  search  type  request  at  Q-2. 

The  disc  program  search  type  request  at  Q-2. 

The  mean  service  time  for  the  R  requests. 

The  mean  service  time  for  the  R  or  R ,  or  R 

r  as  p  s 

reques t  s . 
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Hie  percentage  of  nit  requests  entering  Q-  1  that 
are  ith  priority  requests. 

The  mean  utility  of  the  computer. 

The  mean  utility  of  the  disc. 

The  mean  number  of  computer  service  units  required 
by  a  new  request. 

Frequency  of  computer  service  requests  generated 
by  a  single  new  request. 

Frequency  of  disc  service  requests  generated  by  a 
single  new  request. 

Mean  length  of  the  M  new”  request. 

Mean  length  of  the  repeat  request. 

The  probability  that  the  first  disc  search  is  a 
data  disc  search. 

The  probability  that  the  ith  disc  search  is  a 
data  type. 

The  probability  that  the  ith  disc  search  is  a 
program  type. 

The  proportion  of  n  disc  search  requests  that  is 
the  B d  ty pe . 

Mean  length  of  n  consecutive  “repeat”  requests. 

Mean  length  of  a  “computer”  request. 

Total  number  of  words  of  all  requests  generated 
by  a  weighted  “new”  request. 

The  weighted  mean  length  per  request  at  Q-  1 . 

The  weighted  mean  length  per  request  at  Q-2. 

Total  number  of  words  of  ail  disc  service  requests 
generated  by  a  new  request. 

The  probability  that  a  request,  upon  its  arrival 
at  Q-  1  ,  finds  that  the  request  that  is  being 
serviced  by  the  computer  has  a  lower  or  equal 
priority. 

The  mean  response  time  of  the  new  requests. 


A.  The  Utility  of  the  Computer 

The  mean  utility  of  the  computer,  Pc ,  is  defined  as  the  percentage 
of  time  that  the  computer  spent  in  servicing  requests.  The  time  that 
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the  computer  spent  in  servicing  the  requests  consists  of  two  parts: 

(a)  the  time  spent  in  the  analysis  of  requests,  and  (b)  the  housekeeping 
time  when  an  interruption  occurs.  These  two  parts  of  computer  utility 
are  represented  by  p1  and  p" ,  respectively. 

The  ^i;o  requests  require  only  one  unit  of  computer  service  for  all 
values  of  i  and  j .  The  R i . n  request,  for  n  >  0,  requires  n  +  2  units  of 
computer  service  if  at  the  time  the  request  arrives  at  Q- 1  the  computer 
is  not  busy,  or  if  at  the  time  the  new  request  arrives  at  Q- 1  the  com¬ 
puter  is  busy  and  its  priority  is  higher  than  that  of  the  request  that  is 
being  serviced  by  the  computer.  Thus  the  proportion  of  Rljn  that  requires 

n  +  2  units  of  computer  service  is  (1  -  p  c)  +  p cZ  *  1  -  ,0^(1  ~  Z).  The 

rest  of  the  Rljn,  Pc ( 1  “  2)  percent,  requires  n  +  3  units  of  computer 
service. 

The  probability,  Z,  the  sum  of:  (a)  the  probability  that  the  new 

request  is  the  i  =  1  type  and  the  request  that  is  being  serviced  is 
either  the  i  *  2  or  3  type,  and  (b)  the  probability  that  the  new  request 

j  =  2  type  and  the  request  that  is  being  serviced  is  the  j  *  3  type.  Hence, 

z  =  YX(Y2  +  y3)  ♦  y2T3  .  (a-d 

Since  X  percent  of  the  new  requests  are  the  ^l;o  type  and 
(1  -  X)  percent  are  the  ^ijn  tyPe»  the  mean  number  of  computer  service 
units,  jV,  required  by  a  new  request  is 

N  -  X  +  (1  -  A){[l  -  pc(  1  -  Z)](n  +  2)  +  yOc(l  -  Z)(n  +  3)} 

-  1  ♦  (1  -  X)  [n  ♦  1  +  pe (1  -  Z)]  .  (A- 2) 

The  mean  service  time  for  each  unit  of  computer  service  is  hc , 
hence  the  mean  service  time  for  each  new  request  is  hcN.  For  an  arrival 

rate  of  X,  the  utility  of  the  computer  for  analyzing  the  request  is 

P‘e  *  XheN 

*  Mjc{ 1  ♦  (1  -  X) [n  +  1  +  pe(i  -  Z)]}  .  ( A- 3 ) 
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Each  new  request  will  cause  an  interruption  of  the  computer  service 
if  upon  its  arrival  at  Q- 1  the  computer  is  servicing  another  request. 

The  probability  that  a  new  request  finds  the  computer  husy  is  p  c  .  For  A. 
new  request  arrival  rate  and  h h  mean  housekeeping  time  per  interruption, 
the  time  spent  by  the  computer  in  housekeeping  processing  of  new  request 
is  '  Each  new  request  will  on  the  average  generate  N  -  1  disc 

searches,  where  N  is  given  in  Eq.  (A -2).  A  disc  search  return  request 
may  interrupt  the  computer  service  if  the  priority  of  the  return  request 
is  higher  than  that  of  the  request  that  is  being  serviced.  The  proba¬ 
bility  that  a  return  request  finds  the  computer  busy  is  again  pc ,  and  the 
probability  that  the  return  request  is  higher  than  that  of  the  request 
that  is  being  serviced  is  Z,  as  given  in  Eq .  (A-l).  Thus,  for  a  new  re¬ 

quest  arrival  rate  of  A,  the  return  request  interruption  frequency  is 
K Pc(N  -  1)Z.  The  time  spent  by  the  computer  in  housekeeping  processing 
of  the  return  request  is  then  ^hhpc(N  -  1 )Z.  Consequently,  the  total 
housekeeping  load  at  the  computer  is 

p"c  --  Mi hPcll  ♦  (N  -  1)2]  .  ( A- 4 ) 

Combining  Eqs.  (A-3)  and  (A-4)  one  has 


Pc  -  Mcil  ♦  (1  -  A)  [n  ♦  1  ♦  pe(l  -  Z)]} 

+  hhPeil  *  (1  -  X)[n  *  1  ♦  pt(l  -  Z)]Z }  .  ( A- 5 ) 

By  rearranging  terms  in  Eq.  (A-5),  one  obtains  a  quadratic  equation 
in  Pc ,  A p^  +  Bp c  +  C  =  0,  where 

A  -  Khh(l  -  A) ( 1  -  Z)Z 

B  -  \{/ie(l  -  A)  ( 1  -  Z)  *  hh[l  ♦  (1  -  A)(n  ♦  1  )Z]  -  1} 

C  =  khc  [l  ♦  (1  -  A)(n  ♦  1)] 

Therefore , 

-B  ±  vV  -  4 AC 


70 


However,  for  A  <<  B  or  C,  which  is  the  case  under  this  investigation, 
Eq .  (A- 6)  may  be  approximated  by 


B.  The  Utility  of  the  Disc 

There  are  four  types  of  requests  that  require  the  service  of  the 
disc;  these  are  the  store  type  requests  R$,  the  pop-up  request  retrieval 
type  R  the  disc  data  search  type  requests  Rdt,  and  the  disc  program 
search  type  requests  Rp g .  As  far  as  the  utility  of  the  disc  is  concerned, 
the  Rr f  and  Rds,  and  R  p g  requests  may  be  treated  as  a  single  type,  their 
service  time  being  exponentially  distributed  with  mean  service  time  h d> 

The  mean  service  time  for  R  is  h  . 

s  s 

The  R  is  generated  when  the  new  request,  which  requires  disc  search, 
upon  reaching  Q- 1  finds  that  the  computer  is  busy,  and  that  the  priority 
of  the  new  request  is  equal  to  or  lower  than  that  of  the  request  in  the 
computer.  The  percent  of  new  requests  that  require  disc  search  is  1-A, 
the  probability  of  finding  the  computer  busy  is  ,  and  the  probability 
of  finding  the  request  being  serviced  by  the  computer  has  a  higher  or 
equal  priority  is  (1  -  Z  ).  Therefore,  the  probability  of  a  new  request 
generates  R  at  Q-  2  is  (1  -  X)pc(l  -  Z)  .  The  disc  utility  due  to  /?j  per 
new  request  is  (1  -  X)pe( 1  -  Z) h ( .  The  number  of  Rr,  Rdg,  and  R  gen¬ 
erated  by  a  single  new  request  is  N  -  X,  where  N  is  given  by  Eq.  (A-2). 

The  disc  utility  due  to  the  R  r ,  R d $ ,  and  R  $  per  new  request  is  (N  -  X) h  d . 
Therefore  the  total  disc  utility  is,  for  a  new  request  arrival  rate  of  X, 

Pd  =  Ml  -  X){[n  +  1  +  pc{  1  -  Z)]hd  +  pe(  1  -  Z)ht)  .  ( A- 8 ) 

C.  The  Mean  Hesponse  Time 

The  mean  response  time  of  a  new  request  Tr  is  defined  as  the  mean 
elapsed  time  between  the  arrival  of  the  new  request  at  Q- 1  to  the  time 
when  it  leaves  the  computer  after  the  final  disc  search  has  been 
comp  1 e  ted . 

From  the  response  time  point  of  view,  the  new  requests  may  be  clas¬ 
sified  into  three  categories.  The  first  is  the  requests  that  need  no 
disc  search  R^  X  percent  of  the  requests  are  of  this  category.  The 
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second  is  the  R[ }  n  type  of  request  that  finds  the  computer  not  busy  when 
it  arrives  at  Q-l  or  finds  the  computer  busy  but  its  priority  is  higher 
than  that  of  the  request  in  service.  Thus,  the  percent  of  new  requests 
that  are  in  the  second  category  is  (1  -  ,Y)[(1  -  p  )  +  P CZ]  .  The  third 

is  the  R"  n  type  of  request  that  finds  the  computer  busy;  its  priority 
is  lower  than  or  equal  to  that  of  the  request  in  service.  The  proportion 
of  new  requests  that  belong  to  the  third  category  is  thus  (1  -  Z)  . 


The  mean  response  time  of  the  first  category  of  request  Tj  is  the 
mean  waiting  time  plus  the  mean  computer  service  time  hc  for  the  re¬ 
quest.  In  this  case,  the  mean  waiting  time  is  that  of  a  single¬ 
exponential  server  system,  namely 


Therefore, 


P'K 

1  -  P' 


h  1  ♦ 

c 


(A-9) 


The  second  category  of  requests  require  n  +  1  successive  services 
at  the  disc  and  the  computer.  Each  service  at  the  disc  or  the  computer 
consists  of  a  waiting  time  in  Q- 2  or  Q- 1  and  a  service  time  in  the  disc 
or  the  computer.  The  service  time  of  the  disc  is  exponentially  dis¬ 
tributed  with  mean  service  time  hd,  and  the  disc  is  a  single  server. 
Thus,  the  mean  delay  experienced  by  a  request  at  the  disc  is  that  of  a 
s l ng 1 e- exponent i a  1  server,  or  h d/ { 1  -  p d) .  The  mean  delay  experienced 
by  the  request  at  the  computer  is  as  given  in  Eq.  (A-9).  Since  Q-  1  and 
Q-  2  are  unlimited  in  size,  the  Q- 1  and  Q- 2  are  statistically  independent 
of  each  other.  Thus,  the  mean  response  time  for  the  second  category  of 
requests  ?2  is  merely  n  +  1  times  the  sum  of  the  mean  delay  at  the  disc 
and  the  mean  delay  at  the  computer,  or 
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2 


(1 


In  +  1) 


+ 


(A- 10) 


The  third  category  of  requests  requires  n  +  2  successive  services  at 
the  disc  and  the  computer.  in  addition  it  requires  one  additional  service 
at  the  disc,  the  storing  of  the  request  in  the  disc.  The  mean  service 
time  of  this  last  type  of  disc  service  is  h  .  Thus,  the  mean  response 
time  of  the  third  category  of  requests  T 3  is  n  +  2  times  the  sum  of  the 
mean  delay  at  the  disc,  the  mean  delay  at  the  computer  plus  one  mean  wait¬ 
ing  at  Q- 2,  and  one  disc  service  time  whose  mean  is  h  for  storing  the 
request,  or 
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in  *  2) 
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(A- 11) 


The  total  response  time  Tr  is  the  weighted  sum  of  response  time  of 
each  of  the  three  categories,  or 

Tr  =  AT j  ♦  (1  -  A )  [  1  -  pc  *  PcZ]T2  +  (1  -  X)pe{l  -  Z)T3 


Xh 


1  -  ft 

'  c 


♦  (1  -  X)  (n  ♦  1) 


1  “  P. 


Pa 


♦  PM  -  Z) 


h.  +  h 

d  i 


1  -  Pr  1  -  P. 


(A- 12) 


D.  The  Mean  Response  T  i  m e  of  ith  Priority  Request 

In  this  service  system  the  incoming  requests  are  categorized  into 
three  priority  classes.  All  incoming  requests  require  at  least  one  initial 
service  by  the  computer.  1  -  A  percent  of  all  incoming  requests  require 
n  +  1  subsequent  services  by  the  disc  and  the  computer.  Of  the  1  -  X  per¬ 
cent  requests,  those  requests  that  find  the  computer  serving  a  higher  or 
equal  priority  request  when  they  arrive  at  the  computer  require  an  addi¬ 
tional  service  by  the  disc  and  the  computer.  Since  the  probability  of 
finding  the  computer  busy  is  Pc  and  t he  p robab l 1 i t y  of  an  i  or  higher 

priority  request  is  being  serviced  is  £  Y  ,  the  probability  of  an  ith 

j  -  \  j  i 

priority  request  requiring  an  additional  service  is  (1  -  X)p  Z  Y  .  The 

c  j  =  l  ^ 

number  of  subsequent  computer  services  required  by  an  ith  priority  request 
i  s  ( 1  -  A' )  ( n  +  1  +  p  2  Y  ) . 

e  j  =  1  ; 

Let  i  =  1,  2,  3,  be  the  first  request  for  the  service  of  the  com¬ 

puter  by  the  ith  priority  requests  and  R  t ,  i  =  4,  5,  6,  be  the  subsequent 
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request  For  i  ho  service  of  the  disc  and  (  he  computer  by  the  ith  priority 
requests.  bet  /*  be  the  loading  oF  the  computer  due  to  the  ft..  I  f  A.  is 
the  mean  arrival  rate  oF  the  incoming  requests  and  )  t  is  the  proportion 


the  requests  that  are 

of  the 

ith  priority, 

t.  Ii  e  loading  d  u  e 

to  R  is 

i 

P  =  Kh  Y 

l  C  l 

i  = 

1.  2,  3 

(A-  13) 

> 

II 

s- 

1 

fn  +  1  ♦ 

'•  %  r) 

i  =  4,  5,  6 

(A- 14) 

In  this  analysis  both  the  computer  and  the  disc  are  assumed  to  be 
exponential  servers.  The  queue  discipline  of  the  computer  is  the  in¬ 
terrupt  priority  queue  type,  a  Ri  can  interrupt  the  service  of  R  ^  if 
1  j.  The  queue  discipline  of  the  disc  is  the  ordered  queue  type.  The 
mean  time  spent  by  a  R  .  request  in  waiting  for  and  service  by  the  computer 
rs  given  as2 


T' 

i 


(i  -  pi_l)( 1  -  p, i 


( A- 1 5  ) 


where  h  rs  the  mean  service  time  of  the  computer  and  P  =  £  p.. 

c  1  k  =  l  1 

The  time  spent  by  a  request  in  waiting  for  and  service  by  the  disc 
is  not  distinguished  by  the  priority  type  (the  queuing  discipline  at  the 
disc  is  the  ordered  queue  type),  but  by  the  type  of  disc  service  required 
From  the  service  time  point  of  view  there  are  two  basic  types  of  disc 
service  requests,  the  disc  search  type  and  the  disc  store  type.  The  mean 
time  spent  by  either  of  the  two  types  of  requests  in  waiting  for  the 
service  by  the  disc  is  ( )  /  (  1  -  pd  )  ,  which  is  the  mean  waiting  time  of 
a  single  exponential  server  w i t b  mean  loading  of  p d  and  mean  service  time 
of  h d .  Let  the  mean  service  time  of  the  disc  search  request  and  the  disc 
store  request  be  hgf  and  hgp  respectively.  The  mean  time  spent  by  a  disc 
search  request  and  a  disc  store  request  in  waiting  for  and  service  by  the 
d i s  c  are 


T 


d 


h 


s  e 


+ 


hdPd 


( A- 1 6  ) 
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a  nd 


hiPd 

h  + - 

"  i  -  p. 


(A- 17 ) 


respectively. 


An  ith  priority  incoming  request  require  an  initial  service  by  the  com¬ 
puter,  the  mean  response  time  of  this  service  is  T'.,  1  -  X  percent  of  the 

ith  priority  incoming  requests  require  n  +  1  subsequent  services  by  the 
computer  and  n  +  1  disc  search  type  service  by  the  disc.  The  mean  response 

time  for  the  total  of  n  +  1  services  by  the  computer  and  the  disc  is 

i 

(n  +  1)(T*  ,  +  T,).  In  addition  (1  -  X  )/0  £  Y.  percent  o  f  t  he  1 1  h  p  r  i  o  r  i  t  y 

1  J  a  C  J  c  i  J 

request  requires  a  subsequent  service  by  the  computer  and  a  disc  store  type 
service  by  the  disc.  Thus,  the  grand  total  response  time  of  an  ith  priority 
request  is 


r 


(i  -  x)  [(n  +  i ) ( r ;  +  3  ♦ 


Td) 


Pc 


1 

/- 1 


Y  ( T ' 

; '  i  + 


Tp)] 


(A- 18 ) 


E.  The  Mean  Storage  Size 

The  mean  storage  size  is  defined  as  the  average  number  of  words 
stored  in  Q- 1  and  Q- 2.  Since  the  computer  and  the  disc  are  the  exponen¬ 
tial  type  service  system  the  mean  number  of  requests  inQ-lisL^j*  1/1  -  P  c 
and  in  Q- 2  is  =  1/1 p d.  The  mean  number  of  words  in  Q- 1  and  in  Q- 2  are 


Lq  1  =  1 

•  h  * 

i,/( i  - 

Pc') 

( A- 19 ) 

/  s  / 

LQ 2  LR 2 

'  ~2  = 

V<  i  - 

p rf) 

( A- 20 ) 

respectively,  where  Pc  and  p 

are  given 

by  Eqs. 

( A- 7 )  and  (A-8). 

L ,  and 

L2  are  the  mean  length  of  a  request  in  Q- 1  and  Q- 2  respectively. 

Before  proceeding  to  find  the  L ^  and  L 2,  the  relative  frequency  of 
occurrence  and  the  mean  length  of  various  types  of  request  are  to  be  ex¬ 
amined  from  the  storage  point  of  view.  New  requests  can  be  classified 
into  three  categories,  namely  R . .  Q ,  R[jn  and  and  their  relative 

frequency  of  occurrence  are  X,  (1  -  [l  -  pc  ( 1  -  Z )]  and  (  1  -  X)pc  ( 1  -  Z) 
respectively.  Each  of  these  three  “new"  requests  generate  other  types 
of  requests  as  they  are  processed  by  the  computer  and  the  disc.  The  R.^ Q 
request  generates  a  “new"  request,  R n ,  at  Q- 1  only.  The  R[  ■ n  request 
generates  a  Rn  at  Q- 1  and  at  the  same  time  generates  n  +  1  disc  search 
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requests ,  R  ,  at  Q-2.  Each  of  the  first  n  It  s  requests  hi  term  generates 
a  11  repeat*'  request,  Rp,  for  Q  1.  The  last  R  s  request  generates  a 

“computer"  request,  It  c  ,  for  Q-  1  The  fl  “  request  generates  a  Rn  at  Q-  1 

and  at  the  same*  time  generates  a  disc  store  request,  Rd>  for  Q-  2  and  a 
pop-up  list  request,  Rp ,  for  Q- 1 .  When  the  is  processed  by  the  com¬ 

puter  a  pop-up  list  retrieve  request,  R pr,  is  generated  for  Q- 2,  which 
results  in  the  retrieval  of  the  stored  ft'  from  the  disc  and  is  placed 

into  Q-  1 .  The  retrieved  R'n  is  processed  as  the  Rn  generated  by  the  fl  [• 

request.  Table  A- 1  is  a  summary  of  the  number  of  each  type  of  request 
generated  bv  the  R  ftf  R *  and  R"  request. 

°  i  j  0  *  i  j  n  i  ]  n  * 


Table  A- I 

NUMBER  Of  GENERATED  REQUESTS 


ORIGINAL  REQUEST 

GENERATED  REQUEST  TYPE 

Type 

Relative  Frequency 
of  Occurrence 

Q- 1 

<?-  2 

R 

n 

rx 

R 

P 

R 

r 

R 

c 

R 

Pr 

s 

R  n 

ij  0 

X 

1 

/?' 

ijn 

(1  -  A)  [l  -  pe(  1  -  Z)] 

1 

n 

1 

n  +  1 

R" 

ijn 

(1  -  X)pe ( 1  -  Z) 

1 

1 

1 

n 

1 

1 

1 

n  +  1 

Let  F  ,  F'  ,  F  ,  F  ,  F  ,  F,,  F  ,  and  F  be  the  number  of  the  R '  ,  R 

n  n  p  r  c  d  pr  s  n  p 

Rr,  Rc  j  Rd,  /?  ,  an(l  Rs  requests  generated  by  an  original  request  respec 

tively.  From  Table  A- I  it  is  seen  that 


Fn  -  X  ♦  (1  -  X)[l  -  p( 1  -  Z)]  ♦  (1  -  A)pc(l  -  Z)  =  1 

and  similarly 

K  -  ^  -  r,r  -  (i  -  *>pcu  -  2) 

Ff  -  n ( 1  -  A) 

F  =  1  -  A 

C 

Ft  =  (n  ♦  1)(1  -  A)  .  (A- 21 ) 

Let  £  j  and  £ 2  be  the  number  of  requests,  regardless  of  type  gener¬ 
ated  by  an  original  request  at  Q- 1  and  Q-2  respectively,  then 

fi  =  fn  +  +  Fp  +  rr  ♦  FC  =  1  +  (1  -  A )  [  n  ♦  1  +  2p c ( 1  -  Z ) ]  ,  (A- 22) 
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(1  -  X)ln  ♦  1  ♦  2pe{\  -  Z)] 


( A- 2  3 ) 


a  nd 


F,  +  F  +  F 

d  p  r  s 


Let  A/  ,  M'  ,  M  ,  A/  ,  M  ,  A/.,  and  M  be  the  mean  length  of  the  ft  . 

n  n  p  r'c  a  i  °  n 

ft',  ft  flr,  ftf,  ftrf,  and  ft  requests  respectively.  The  mean  length  of 
ftpr  is  the  same  as  that  of  fip.  The  weighted  mean  length  per  request  at 
Q-  1  is 


(F  M  ♦  F#A/'  ♦  F  A/  +  F  M  +  F  Af  ) 

'  n  n  nn  P  P  r  r  cc 


( A- 24 ) 


and  at  Q- 2  i 


+  fpA  + 


( A- 25 ) 


Only  the  Afn  and  A/p  are  the  basic  system  parameters.  The  other  mean 
length  of  requests  are  to  be  developed  in  term  of  the  basic  system  pa¬ 
rameters,  namely,  A/  ,  A/  ,  V  ,  A/,  ,  and  A/ 

The  expressions  for  Md  and  M*n  are  to  be  developed  first.  When  a 
R'lJn  is  processed  by  the  computer  the  original  request  plus  the  disc 
search  address  for  each  of  the  n  ♦  1  disc  searches,  plus  the  disc  store 
address  are  forwarded  to  Q- 2  as  a  Hd .  The  mean  length  of  each  disc  ad¬ 
dress  is  A/  ,  therefore  the  mean  length  of  Hd  is 

A td  =  Mn  ♦  (n  +  2)A/p  .  (A-  26) 

When  the  ft^  is  processed  by  the  disc  the  request  minus  the  disc  store 
address  is  stored  into  the  disc.  This  request  when  retrieved  by  a  Rpr 
becomes  a  ft '  ,  thus  the  mean  length  of  R '  is 

n  °  n 


K  *  Md  -  Mp  =  Mn  ♦  (n  +  l),Wp  .  (A- 27) 


A/c  is  the  mean  length  of  the  display  message  or  R  retrieved  by  the 
last  R t  from  the  disc.  There  are  three  types  of  display  messages,  their 
proportion  and  mean  length  are  Vi.  and  l  t  j  *  lt  2,  3,  respectively. 
Hence,  the  mean  length  of  ft  is 


M 


1  w  T 

j  =  i  11 


(A- 28 ) 
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To  dove ! op  the  expression  for  Mr  and  M s  it  rs  necessary  to  review 
the  process  of  handling  the  /?  and  R  requests  by  the  disc  and  the  com¬ 
puter.  The  first  disc  search  request,  fl  j ,  has  a  search  length  of 
L  l  =  (n  +  l )  \!p .  When  /i  j  is  processed  by  the  disc  a  message  is  retrieved 
from  the  disc.  This  message  and  the  Rrl  (its  search  length  reduced  by  one 
M  unit)  are  placed  into  Q-  1  to  form  the  first  "repeat”  request,  ^p\’ 

Thus,  the  length  of  R  ,  is  /-  ,  =  nM  +  M  where  M  is  the  mean  length  of 
the  retrieved  message.  The  R f  2  is  processed  by  the  computer  and  is  placed 
into  ()-2  to  become  R  g  0  .  From  reques  t  -  1  eng  t  h  point  of  view  L$2  -  LrJ. 

When  the  R  is  processed  by  the  disc  the  message  retrieved,  plus  its 
search  length  reduced  by  one  M  unit,  are  placed  into  Q-  1  to  become  the 
R  p  2  •  The  length  of  Rp2  is  then  lf2  =  ( n  -  1  )Mp  +  .  The  2  is  processed 

by  the  computer  and  is  placed  into  Q-2  to  become  R 3.  This  process  is  re¬ 
peated  n  times.  In  general,  the  length  of  Rpr  is  Lrk  =  ( n  +  1  -  k)  M  p  +  > 

and  Lrk  -  +  for  k  ^  1.  Thus,  the  mean  length  of  a  R p  is 


n  n 

»,  ■  1  Z_  1r.  *  1  Z  l"  * 

M  *  n  1 


1 

y 

n 

h  =  1 

mean  1 

1 

n 

♦  1 

re 

a  re 

1  -  fe):W  +  M  ] 

P  * 


n  +  1 

-  M  +  A/  (A- 29) 

2  p  * 


k  ~  l 


i  +  l 


k  =  l 


_ 1_  ( n  ♦  l)in  +  2) 

n  +  1  2 


jW  +  n  M 

P  a 


(A- 30) 


program  disc  search  type,  /?  .  The  mean  message  length  retrieved  by  R  d  g 

is  A/^s,  and  that  retrieved  by  R  is  M  g  .  Let  F^  be  the  proportion  of  R 
that  are  the  Rdg  types,  thus  M$  =  ^n^ds  +  (1  ”  Substitute  the 

expression  for  \l  s  into  Eq .  (A-29)  and  (A-30)  one  has 


n  ♦  1 

-  \1  +  F  M,  ♦  ( 1  -  F  )A/ 

2  p  n  d  s  rx  p  $ 


an  (1 


V 


1  f(n  + 

n  1  1  I 


1 ) (n  +  2) 


«  -  n[K  AC.  ♦  (1  -  V)M 


d 


(A-31 ) 


(A- 32) 
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The  expression  for  Vn  is  a  function  of  the  rules  that  determine  the 

I  * 

type  of  disc  search  for  the  kt\\  disc  search.  The  rules  are  as  follows: 


(1)  the  probability  of  Rgi  being  a  R d g  is  U  and  being 

(2) 


a  R  is  1  -  U, 

p  * 

i  f  ft  .  i  s  a  R  then  R  ,  L  .  ,  v  must  be  a  R,  ,  but  if 

»  k  p  S  S  (  k  T  l  )  d  s 

R  L  is  a  Rdg  then  Rs(k+\)  can  be  a  Rdg  with  a 


probability  of  U  and  be  a  Rpt> 

(3)  the  Rtn  is  always  Rps- 

Let  Pdk  and  Ppk  be  the  probability  that  the  feth  disc  search  is  a  Rdg  and 
R  p  a  respectively.  Since  a  disc  search  can  only  be  a  Rjat  or  a  Rpt> 
therefore  Pdk  3  1  -  Ppk-  According  to  Ru 1 e -  1  Pd  j  *  U  and  Ppi  *  1  -  U. 
According  to  Rule-2  Pd(k  +  j ) *  ^ d  k^  ~  ^  ^ dk ^  “  tO  •  SinceP^  *  1  -  P dk, 

there  fore 


P  (k  +  1) 

p 


Pdk(l  -  u ) 


( A- 33 ) 


Equation  (A- 27)  is  a  recurrence  equation,  once  P d  \  is  given  Pdk  and  Ppk 
can  be  uniquely  determined.  The  proportion  of  the  n  disc  search  that  are 
the  Rdg  type  is 


n  -  1 


k  «  1 


Pdk 


( A- 34 ) 


It  is  to  be  noted  that  the  summation  is  over  n  -  1  since  by  Rule  3  the 
R  is  always  a  R 

s  n  y  p  t 

Substituting  Eqs .  (A-21),  (A-22),  (A-23),  (A-26),  (A-27),  (A-31), 
and  (A- 32 )  into  Eqs.  (A-24)  and  (A- 25),  and  then  substituting  Eqs.  (A-24) 
and  (A-25)  into  Eqs.  (A-19)  and  (A-20),  we  get 


K  +  (1  -*>(" H,  +  VMit  +  (1  -  Vn)Mp, ]  ♦  Mc  ♦  LWn  ♦  (n  ♦  2)Mp)pt{l  -Z) 


"<?l 


(1  -  p c ) {l  ♦  (1  -  JOln  ♦  1  ♦  2pe(  1  -  Z)]> 


(A- 35 ) 


The*e  rule*  ere  in  accordance  with  that  atated  in  tha  program  liating  of  the  "473-L  Simulation** 
Monthly  Statua  Report.  The  main  text  of  the  above  cited  report  aimply  atated  that  tha  probability  of 
a  baing  Rpt  la  t /.  In  order  that  the  analyaia  conform*  with  the  actual  aimulation  modal,  the  rulea 
aa  atated  in  the  program  liating  ara  uaed* 
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( A- 36 ) 


a  ihJ 


-  —  *  ~  M  *n[VM,  +(l-V  ),»/  j  +  lM *(n  +  3)A/  K(l-Z) 

p  n  a  *  n  pi  n  pc 


'Q2 


(1  -f>d)ln  +  1  ♦  ^c(l  -  Z)J 


where  .1/  and  l/n  are  given  in  Hqs.  (A- 28)  and  (A- 34). 

In  the  simulation  program  listing  of  the  11  473- L  Simulation”  Report  1, 
it  was  stated  the  mean  length  of  due  to  a  /f*  request  only,  is  the 

same  as  the  mean  length  of  R ”,  namely,  LJ  j  =  A/n  +  (n  +  1 )  Af  .  We  note 
also  this  mean  length  is  one  M  unit  longer  than  that  used  in  deriving 
Fqs.  (A- 35 )  and  (A-36).  Furthermore,  this  one  unit  of  M n  is  not  dropped 
from  subsequent  Rp ,  Rs,  and  Rc  requests.  For  example,  the  mean  length  of 
R  .  is  L  .  =  (n  +  1  -  k)M  +  M  +  M  ,  and  the  mean  length  of  R  is 
Mc  =  5  W .  I  +  A/^ .  Using  the  new  storage  rule,  which  shall  be  called 
Hule  £  to  distinguish  it  from  the  Rule  1  used  in  deriving  Eqs.  (A- 35 )  and 
(A-36),  the  expressions  for  the  mean  number  of  words  in  Q- 1  and  Q- 2  are 


~Q  1 


+  -y  Hp  +  I'Ai  +  d-rn),Wp,]  *MC  ♦  (n  +  2)(Mn  +  Mp )pc ( 1  - Z)| 

( 1  -Pc ) {1  +  (1  -  A)  [n  +  1  +  2pc(l  -  Z)]} 


( A- 37 ) 


a  n  d 


(n  +  l)(n  +  2) 


Vnll,A,  +  +  <#„  *^)(n*3)pe(l  -  Z) 

(1  -  p.)[n  ~  1  ♦  2pJ[  -  Z)] 


( A- 38 ) 


It  IS 
( A- 38 )  are 
results  of 


believed  Rule  2  is  not  logical;  nevertheless,  Fqs.  (A -37) 
presented  here  so  that  they  may  be  used  to  check  with  the 
simulation. 


and 
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APPENDIX  B 


CONFIDENCE  INTERVALS  FROM  A  SIMPLE  SAMPLE 

1  INTRODUCTION 

Assuming  sample  sizes  from  the  simulation  are  large  enough  to  make 
average  output  approximately  normally  distributed,  symmetric  confidence 
intervals  are  desirable  measures  of  variation  in  the  estimates.  This 
section  is  primarily  concerned  with  obtaining  the  confidence  interval. 

A.  Confidence  Intervals 

In  determining  confidence  intervals,  assumptions  are: 

(1)  The  transition  probabilities  are  ergodic.  This  is 
necessary  because  only  one  simulation  run  will  be 
made  from  a  fixed  initial  state. 

(2)  The  initial  state  can  be  chosen  so  that  every  sub¬ 
set  of  states  having  a  substantial  stationary  prob¬ 
ability  will  be  readied  with  high  probability  during 
feasible  sampling  periods  n.  This  insures  that 
within  the  sampling  period  n,  the  stationary  distribu¬ 
tion  of  states  is  approximated  without  leaving  out 
important  subsets  of  states.  Work  is  needed  to  make 
this  requirement  more  precise  in  terms  of  the  transi¬ 
tion  probabilities  of  the  underlying  model. 

(3)  The  autocorrelation  of  the  output  variable  is  geometric 
or  of  the  “decay''  type,  i.e.,  if 

E{X2t)  =  cr2  +  fi2  and  E(X  (X  f+  =  L*J2  ♦  fi2 

then 

E(*,Xt.K)  -  Pkcx2  ♦  fj.2  (B-l) 

This  is  more  flexible  than  assuming  no  autocorrelation  and  it  agrees 
with  observed  output  in  several  queuing  models.  The  validity  of  it  can 
be  tested  by  estimating  £(XtA<+fc)  for  k  =  2,  3,  ...  .  This  assumption 

will  not  hold  when  periodic  effects  occur. 
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(1)  At  least  100  values  of  X f  are  recommended.  Compensation 
for  small  sample  sizes  is  not  included  in  the  procedures. 
Normally  n  <  100  anyway. 

The  output  time  series  Af,  1  ^  f  £  n ,  i s  treated  apart  from  other 
output  time  series.  It  is  not  necessarily  the  raw  output,  but  usually 
is  an  average  of  L  successive  output  values  yfl,  y 1 2  >  •••>  y  tL : 


L 


L 


*  y 

j  =  i 


‘  i 


(B-2) 


This  reduces  computer  time  for  statistical  computation  which  is  propor¬ 
tional  to  1  /L.  \alues  of  L  ~  5  or  10  are  recommended.  As  an  upper 

restriction  on  L,  it  should  not  be  desirable  to  stop  simulation  before 
100  groups  of  L  events  have  occurred. 


Some  definitions  are  needed. 

2  (  y \  rr2  | 


o 2  (A)  and  cr2  (X)  are  the  variances  of  X  f  and  A. 


S2  (A)  and  S2  ( A)  are  sample  variances  of  Af  and  A. 


A\ 

n  1 


—  £  X  X  k 

L  L  t  •  An  2 


k  t~k 


i  n  -  k 
1  V  V 

n-k  «=1  A‘ 


( B-  3 ) 


n-k 

1 

t  =  i 


-  *„*,)<*.♦*  -  xnk2) 


(B-3a) 


p\ 


(n  -  k)S  2(X) 


is  the  sample  order  autocorrelation  .  ( B-  3b) 


n  -  2k  '  1  +  Pk 
(n-k)2  ( 1  -  pk  )2 


for  K  =  1,  2,  3, 


(  B  -  3  c  ) 


Then  the  confidence  interval  about  X  is  X  i  aS(X)  where  a  is  deter¬ 
mined  by  the  level  of  confidence  level  and 


S2(X) 


n  -  l 


2  I 

*  =  i 


(n  -  k)p* 


< 


s2(x)  n  +  pe 

n  1  -  PE 


Pe  --  Pi  +  H  ~P^D i- 


(B  4) 
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Foi  large  n,  this  is  almost  an  equality. 

Before  obtaining  a  confidence  interval,  one  needs  a  "reliable" 
estimate  of  /  .  It  is  found  that 


P  £  =  P  J  +  (1  —  ^  1  )  ^  1 


(B-5) 


is  almost  unbiased,  especially  when  D ^  is  small.  Criteria  considered 
for  a  good  estimate  of  p  are 


eC  *  p,\ 

1  -  p 

(\  ♦ 
s  !  I 

a . 

V  -  pj 

1  -  p 

[<  - 

P, 


h  >  0 


( B-6a ) 


1  ♦  P 


1  -  P 


h ' 


1  ♦  Pt 


0.5, 


h '  >  1 


(  B- 6  b ) 


If  h  is  small,  bias  is  small.  If  h'  is  small,  variation  is  small. 

When  £)j  <0.1  and  n  >  100,  it  is  probably  true  that  h  <  0.1  and 
h‘  <  1.5,  although  no  proof  of  this  has  been  found. 

To  see  why  D l  is  involved,  consider 


E{Ck) 

(n  -  k)o2{X) 


(B-7) 


and  particularly  consider  k  =  1. 


flSi 

(n  -  k)o2(X) 


n__2k  /l  ♦  p\ 
(n  -  k)2  V1  -  P/ 


terms  in 


+  terms  in 


(n  -  *)2  j 


k  =  1,2,3,... 


( B-  8 ) 


For  n  >  100,  and  k  small,  the  last  two  terms  are  insignificant.  The  bias 
is  in  the  second  term  which  can  be  approximated  by  (1  -  Pl)Dl  However, 
since  the  bias  does  depend  on  p  which  we  are  estimating,  a  stable  solu¬ 
tion  of  this  equation  occurs  only  for  small  values  of  the  second  term 
relative  to  1  -  p.  D j  is  a  sample  estimate  of  the  second  term  over  1-p 
In  this  way,  D ^  and  h  are  related. 


I  * 


83 


Further,  for  A'  normally  distributed  an  approximate  0  95  probability 
bound  on  p  was  obtained  for  cases  where  D ^  0.1.  It  is  shown  in  the  next 

section  that  h '  1.5  except  when  n  <  100.  This  partially  substantiates 

the  relationship  between  Dl  and  h'  . 

B.  Variance  of  the  Correlation 

The  stopping  rule  depends  upon  finding  a  relatively  unbiased  estimate 
of  p,  the  correlation.  Whether  this  provides  an  estimate  of  [(1  +  P)/(l  -P)] 
with  small  variance  is  the  critical  question.  This  section  shows  that  for 
D  =  0.1,  and  f  >  1/2,  assuming  normally  distributed  output, 


<  0.05 


with  h  = 1.5 


(B  9) 


This  would  seem  to  provide  a  reliable  estimate  of  [(1  +  p)/(l  -  p)]. 
Note  that  h  is  a  function  of  the  stopping  rule  and  D  in  particular. 
Other  values  of  D  could  be  used,  but  only  Z)  =  0.  1  is  investigated  here. 

First,  it  can  be  reasonably  assumed  that 


X 

l 


L 


L 
1 
j  -  i 


i  ) 


is  normally  distributed,  for  large  values  of  L.  For  L  =  5  or  10,  this 
is  frequently  not  a  very  good  assumption,  but  this  seems  to  be  the  only 
case  for  which  the  distribution  of  correlation  has  been  analyzed. 

Then  the  key  to  this  development  is  a  result  of  Fisher  that  for 
normal  variates  X  and  Y,  Z  =  1/2  log^  [(1  +  Pjy,)/(1  -  Pjy)]  is  normally 
distributed  with  variance  cr^  =  \/(n-  3),  where  n  is  the  sample  size.  In 
this  case  each  pair  (A  ),  t  =  1,2,  ..  n  is  independent  of  the  others. 

In  serial  correlation,  F  t  =  X t  +  ^ .  The  estimate  of  serial  correlation 

is  based  on  correlated  observations,  but  with  approximately  the  same 
distribution  as  Z  with  cr^  =1  (n  -  3)[l  +  pz)/(l  -  P2)]. 

An  estimate  of  p^,  the  correlation  between  X  (X  (  l  and  is 

needed.  It  turns  out  that  pt  <  p£.  Derivation. 

Under  the  geometric  assumption  used  to  derive  p£ ,  Kt  is  given  by 

-  ^,.i  +  e. 
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wli  e  re 


and 


S  „(o.  —  ;).  .11  - 

A  t  "v  n  (  0  ,  cr2  )  a  1 1  t 

=  0,  K  /  0 

-  0  *  /  0 


assuming  mean  0  is  easiest  but  completely  general.  In  general, 
that  £*  ( A  ^  )  =  /fcr4 .  (K  -  3  for  the  normal  distribution.  )  Then 


£Lv,.v,.,: 


=  A/’V  +  ( 1  _  p2  )a4 


+  e,)2(/°2^,.1  +  P6, 


+  6 


(  +  1 


>*,J 


=  E[(p'x\_x  ♦ 

•  (p2*,  ,  ♦  pc,  ♦  6,m)J 

=  Ap4cr$  t  2p2(l  -  p2)cr4  ♦  Jo2(l  -  fj2  )ct4 

=  pV4  [/>2A  +  3(1  -  P2)] 

C  =  -  {ALV,^,^]}2 

=  p2a*(p2h.  *  3)  ( 1  -  p2)  -  P2o-4 

(K  -  3)P2  ♦  2  ' 

(K  -  Dp2  +  1 . 


(B- 10) 

assume 


(B-ll) 
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Win'll  is  f. '  z 


(A  -  3 )P2  +  2  1 

-  <  - 

(A  -  2)p2  *  1  p 
(A  -3)/ 3  _  (A  -  2)P2  +  2p  -  1  <  0 
{f  - ] )  ( A  -  3 ) p2  ♦  1  -  P  <  0,  so  A'  •  ( B-  12) 


For  the  normal  distribution,  A  3  satisfies  this  inequality,  so  that 
estimating  P  z  by  p ^  will  be  conservative.  Note  that  for  p  >  4,  P2  is 
close  to  p ,  and  for  small  pt  we  have  pz  =  p2  .  (For  other  common  dis¬ 
tributions  A  >  3,  i.e.,  for  exponential  K  -  9.) 


Now  given  an  estimate  p the  variance  of  z  is  estimated  by 
(1  +  p£)/[(ri  -  3)(1  -  p^-)]  and  the  upper  critical  point  of  z  for  OL  =  0.05 
i  s 


( B- 13) 


This  can  be  inverted  to  yield  the  upper  critical  value  of  p  called  p. 


P 


(B- 14) 


For  p£  =  0.2,  0.4,  0.6,  0.8,  0.9,  0.95,  these  values  are  given  in 
Table  B- I ,  from  which  h  1.5  is  derived  empirically. 

For  small  p,  a  good  estimate  is  guaranteed  by  requiring  n  >  100. 
For  p  =  0,  as  an  example, 


Pr{pE  >  0.2)  <  0.05  if  n  =  100,  ( B- 1 5 ) 

because  this  yields  100  independent  observations  on  p,  and  this  distribu¬ 
tion  has  been  worked  out  for  normal  variates.14 
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Till* l o  H-l 


\ITHO\IMATK  \'\'\'YM  CHITICM  VAIJ'KS  01*  [(I  ♦  p)(l  “  p)]  KUH  a  *  0.05. 
p  {  [  (  I  +  p)  {  I  -  p  )  *  /i  [  (  1  *  P  e  )  /  (  1  -  P  £ )  ]  }  -  0 .05, 

\s.s  N  nr  in  a  I  Variates. 


A 

\\ 

c 

1  +  E 

 <  fn  10 

■»  _  n 

•i  (j 

U  z 

BC  -  1  - 

1  ♦  P 

♦ 

h 

I  -  Pf 

i .  pk 

n  -  3 

-  z 

w  p 

BC  *  1 

1  -  P 

0 .2 

18.7 

0.1 

0.  64 

1.0 

1.75 

- *  0.47 

3.75 

3.0 

2.0 

1  .  5 

0.  1 

2.33 

37 

0.  On 

0.50 

1.65 

2.84 

-  *  0.50 

4.84 

4.0 

1.7 

O.o 

\ 

100 

0.04 

0.4 

1.5 

5 

-  =  0.71 

7 

5.7 

1.4 

0.8 

n 

150 

0.02 

0.28 

1.33 

11 

—  *  0.85 

13 

12 

1.5 

0.0 

10 

1000 

0.01 

0.2 

1 . 22 

22 

—  --  0.92 

23 

24 

1.3 

O.05 

30 

7800 

0.005 

0. 1  X 

1.15 

44 

—  *  0.956 

45 

45 

1.15 

Not*  that  no  upward  adjustment  for  bias  in  P g  is  made.  Also,  no  downward  adjuatment  is 
made  for  the  known  difference  between  Pg  and  P For  P  ^  1/2,  theae  adjustmenta  are 
cancelling  and  negligible. 


For  small  p,  p  ^  =  Pg  and  any  adjustment  should  be  downward,  ao  that  h  1a  overeat imated 
here . 


1  +  Pe 

Note  that  for  p ~  <  0.2,  -  <  1.5  as  desired. 

4  1  -  Pe 

Further,  even  with  p  -  0,  one  would  expect  correlation  between 
\  Y, +  1  and  Af+1A'  t+2*  But  these  are  independent  when  p  *  0,  so  that  100 
independent  observations  are  obtained  for  n  *  100,  p  -  0. 


C.  Variance  ok  pE  = 


1  /  K 


For  K  >  1,  the  estimate  p *  of  pK  has  the  same  form  and  variance  as 
p  j  does  as  an  estimator  of  p.  But  when  the  Klh  root  is  found,  the  vari¬ 
ance  is  magnified  as  K  increases;  1  K  becomes  a  very  poor  estimator 

of  p.  FIxperience  has  shown  that  this  can  be  a  serious  effect.  Thus,  use 

of  /  ^ \  i  ^ 

\  K)  for  >  1  is  not  recommended  for  small  sample  sizes,  despite 

the  risk  incurred  if  the  geometric  correlation  assumption  is  invalid. 
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D.  Experience  with  Confidence  Interval  Procedure 


It  was  once  felt  that  pj,  p2,  •  ••»  P 5  might  be  good  estimators  of  p. 

However,  we  see  that  only  p ^  is  good  for  small  samples.  This  is  because 
p j,  P2,  P 3 1  •••»  p\  have  the  same  variation.  For  example,  (  ^)i,s  '">s 

more  variation  than  p j.  But  for  large  sample  sizes,  the  validity  of  the 
geometric  assumption  can  be  tested  using  the  sample  estimates  of  higher 
correlations. 

When  only  p ^  was  estimated  and  L  -  5,  all  of  the  statistical  compu¬ 
tation  added  about  20%  to  the  basic  running  time  of  a  simulation.  If 
L  =  .0,  the  computation  added  would  be  about  10%  instead  of  20%,  so  an 
appreciable  gain  is  possible  when  large  grouping  is  feasible. 

This  procedure  has  thus  far  agreed  with  theoretical  results  when 
applied  to  both  input  and  output  variables.  In  several  instances  the 
designed  mean  of  exponential  input  was  verified  within  a  0.034%  interval 
about  the  observed  average.  Similarly,  within  10%  of  the  observed  aver¬ 
age  of  output  variables,  theoretical  values  were  found  as  specified  from 
sample  confidence  intervals.  In  some  cases  the  theoretical  value  is 
close  to  the  boundary  of  the  confidence  interval  so  that  the  intervals 
do  not  seem  to  be  grossly  in  error. 
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APPENDIX  C 


CHOPPING  RULES  AND  BIAS  FROM  INITIAL  STATES 


If  the  initial  state  is  not  selected  at  random  from  the  true  distn- 
bution  of  states  when  the  system  is  in  steady-state  the  choice  will  intro¬ 
duce  some  bias  into  the  resultant  average  statistics.  It  is  desirable  to 
determine  and  possibly  eliminate  this  bias. 


Past  suggestions  have  been 

(1)  Ignore  this  effect  and  start  from  zero,  no  elements 
in  system. 

(2)  Use  a  preliminary  run  and  choose  the  initial  state 

at  random  from  the  end  section  of  the  preliminary  run. 

(3)  Drop  the  first  part  of  the  run  from  average  statistics. 

We  are  using  a  stopping  rule  that  detects  high  correlation,  which 
would  accompany  a  large  effect  from  the  initial  state.  Thus,  the  stopping 
rule  tends  to  force  a  large  sample  to  minimize  the  effect  of  the  initiali¬ 
zation.  Experience  also  indicates  that  we  can  ignore  chopping  rules. 


The  steady  state  probability  of  “idle  state”  occurring  is  easily  de¬ 
termined  as  p  =  (/i  ~  X)  /  p  where  X  is  the  arrival  rate  (all  priorities)  and 
p  is  the  service  rate  (all  priorities).  This  holds  true  except  for  the 
case  of  preemptive  priorities  and  s ta t e- dependen t  arrival  and  service  rates. 
In  these  cases,  conservative  estimates  can  be  found. 


Then  for  samples  of  100  and  p  <  0.Q,  more  than  10  occurrences  of  the 
initial  state  would  be  expected  in  the  run.  Little  bias  would  be  expected 
when  starting  from  such  a  common  state.  So  there  is  justification  for 
ignoring  the  effect  of  initialization.  The  use  of  a  preliminary  run  to 
generate  random  states  produces  several  problems: 

(1)  How  long  does  the  preliminary  run  have  to  be?  This  is  essentially 
the  same  as  asking  about  the  size  of  the  initialization  effect.  Just  knowing 
that  a  preliminary  run  will  decrease  the  effect  is  no  more  information  than 
knowing  that  the  system  will  approach  a  steady-state  condition. 
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( 2 )  If  sovcml  initial  states  are  to  be  recorded  from  the  prelimi¬ 
nary  run,  for  future  use,  a  complete  state  description  of  the  system  must 
be  recorded.  On  the  other  hand,  in  truncating  the  first  part  of  a  run, 
the  state  of  the  system  is  naturally  generated  in  transition. 

The  main  disadvantage  in  chopping  off  the  first  part  of  a  run  is  that 
the  chopping  point  must  depend  only  on  statistics  preceding  the  chopping 
point  or  else  retroactive  adjustments  are  needed. 

We  propose  here  a  method  of  retroactive  chopping  which  r equ 1  res  mi  no r 
retroactive  adjustments  and  which  adaptively  determines  the  chopping  point. 

The  notation  of  the  stopping  rule  will  be  used. 

Assume  that  the  stopping  rule  criteria  have  been  met  for  all  statis¬ 
tical  averages,  so  that  the  run  is  about  to  terminate  on  that  basis.  At 
this  point,  a  check  for  initialization  bias  is  made  on  each  statistical 
average.  After  each  sequence  is  chopped,  the  stopping  rule  is  applied 
again.  If  more  observations  are  needed,  the  same  cycle  will  be  repeated 
until  the  stopping  rule  is  satisfied  immediately  after  chopping. 

Now  consider  the  chopping  procedure  for  one  sequence  X j ,  X 2 . X n. 

Here  n  is  the  number  of  observations  when  the  stopping  rule  is  satisfied, 

—  0  2 

and  the  statistics  Xnt  Sy  and  p£  are  determined. 

n 

Consider  an  intermediate  point  m  where  statistics  X  and  S2  are  stored. 

The  sequence  A^+1,  A^+2,  •  has 

nX  -  mX 

—  n  ■ 

x  =  - 

n  -  at 


[("-OS; 


(m  -  1  )S2  -  (n  -  m)X2  +  nX2  -  mX2]  - - - 

*  n  "J  n  -  m  -  1 

(C-l) 


It  is  unnecessary  to  test  whether  the  mean  fi m  <  Pn_m  significantly 
if  X  m  >  X  But  if  Xm  <  Xn_n,  assuming  equal  variances  for  both  se¬ 

quences,  a  test  statistic  t  normally  distributed  can  detect  significant 
effects  of  initialization.  If  starting  from  the  “idle  state”  would  in¬ 
flate  values  of  X,  then  the  test  would  be  made  in  the  other  direction. 
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Let 


(02) 

If  X  n  _  a  tln<J  X  are  normally  distributed  with  equal  variances,  t  has 
the  {-distribution  with  n  -  2  degrees  of  freedom.  But  for  n  >  100,  this 
is  very  close  to  the  normal  distribution.  So  for  t  >  2,  we  can  reject 
the  hypothesis  that  there  is  no  initialization  effect  and  be  wrong  only 
2.5%  of  the  time.  Because  of  correlation,  this  will  be  a  conservative 
test,  that  is,  t  will  tend  to  be  smaller  as  correlation  increases. 


If  t  >  2,  then  the  sequence  can  be  chopped  at  m,  leaving  A^+1, 
•  ••»  X  n  with  statistics 


X 


P  E 


and 


s2 


7 


+  Pe 


~  Pe 


(C-3) 


Note  that  p£t  as  determined  for  the  full  sequence,  has  not  been  changed. 

If  \Xn_m  -  X m\  <  l/10(Sj  ),  the  chopp  ing  may  not  be  worth  the  trouble. 

n 

This  test  should  not  be  made  often,  because  the  probability  of  chop¬ 
ping  wrongly  tends  to  accumulate  as  a  multiple  of  the  number  of  tests  (7), 
i.e. ,  P(t  >2  in  one  of  T  tests|no  initial  effect)  0.025  7.  It  is  rec¬ 
ommended  that  if  n'  is  the  number  of  samples  between  tests  of  the  stop¬ 

ping  rule,  the  points  m  could  be 

/  / 

n  n  .  .  . 

— ,  — ,  n  ,  2 n ,  4  n  ,  8  n . 

such  that  at  these  points,  Xn  and  5^  are  stored. 
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The  chopping  decision  at  point  m  is  made  only  if  n  ^ 


Th  e  great¬ 


est  value  of  m  with  t  >  2  is  the  chopping  point,  and  the  procedure  con¬ 
tinues  with  the  remaining  n  -  m  points.  This  means  not  only  revising  the 
n-point  statistics,  but  also  scrapping  all  stored  values  at  m-points  less 
than  n,  and  picking  up  new  m  points  as  they  occur  in  the  new  ( n  -  m)  se¬ 
quence.  That  is,  if  n 1  =  100,  n  =  500  and  the  sequence  is  chopped  at  50, 
rt  -  m  =  450  and  all  m  points  are  discarded.  The  first  one  picked  up  will 


be  at  n  -  m  =  800,  and  cannot  be  used  unless  n 


exceeds  1600. 
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